diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json index 4245cf6..4c50ae3 100644 --- a/.vscode/c_cpp_properties.json +++ b/.vscode/c_cpp_properties.json @@ -1,18 +1,16 @@ { - "configurations": [ - { - "name": "Linux", - "includePath": [ - "${workspaceFolder}/ext/pybind11/include/**", - "/usr/include/python3.8/**", - "/usr/include/eigen3/**" - ], - "defines": [], - "compilerPath": "/usr/bin/gcc", - "cStandard": "c11", - "cppStandard": "c++11", - "intelliSenseMode": "clang-x64" - } - ], - "version": 4 + "configurations": [ + { + "name": "windows-gcc-x64", + "includePath": [ + "${workspaceFolder}/**" + ], + "defines": [], + "compilerPath": "gcc", + "cStandard": "${default}", + "cppStandard": "${default}", + "intelliSenseMode": "windows-gcc-x64" + } + ], + "version": 4 } \ No newline at end of file diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..09ca61a --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,24 @@ +{ + "version": "0.2.0", + "configurations": [ + { + "name": "C/C++ Runner: Debug Session", + "type": "cppdbg", + "request": "launch", + "args": [], + "stopAtEntry": false, + "externalConsole": true, + "cwd": "e:/AdaPT/tracer/bvh", + "program": "e:/AdaPT/tracer/bvh/build/Debug/outDebug", + "MIMode": "gdb", + "miDebuggerPath": "gdb", + "setupCommands": [ + { + "description": "Enable pretty-printing for gdb", + "text": "-enable-pretty-printing", + "ignoreFailures": true + } + ] + } + ] +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json index 3ac7a7e..db0591b 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,79 +1,136 @@ { - "files.associations": { - "cctype": "cpp", - "clocale": "cpp", - "cmath": "cpp", - "cstdarg": "cpp", - "cstddef": "cpp", - "cstdio": "cpp", - "cstdlib": "cpp", - "cstring": "cpp", - "ctime": "cpp", - "cwchar": "cpp", - "cwctype": "cpp", - "array": "cpp", - "atomic": "cpp", - "strstream": "cpp", - "*.tcc": "cpp", - "bitset": "cpp", - "chrono": "cpp", - "cinttypes": "cpp", - "codecvt": "cpp", - "complex": "cpp", - "condition_variable": "cpp", - "cstdint": "cpp", - "deque": "cpp", - "forward_list": "cpp", - "list": "cpp", - "unordered_map": "cpp", - "unordered_set": "cpp", - "vector": "cpp", - "exception": "cpp", - "algorithm": "cpp", - "filesystem": "cpp", - "functional": "cpp", - "iterator": "cpp", - "map": "cpp", - "memory": "cpp", - "memory_resource": "cpp", - "numeric": "cpp", - "optional": "cpp", - "random": "cpp", - "ratio": "cpp", - "set": "cpp", - "string": "cpp", - "string_view": "cpp", - "system_error": "cpp", - "tuple": "cpp", - "type_traits": "cpp", - "utility": "cpp", - "fstream": "cpp", - "initializer_list": "cpp", - "iomanip": "cpp", - "iosfwd": "cpp", - "iostream": "cpp", - "istream": "cpp", - "limits": "cpp", - "mutex": "cpp", - "new": "cpp", - "ostream": "cpp", - "sstream": "cpp", - "stdexcept": "cpp", - "streambuf": "cpp", - "thread": "cpp", - "cfenv": "cpp", - "typeindex": "cpp", - "typeinfo": "cpp", - "valarray": "cpp", - "variant": "cpp", - "bit": "cpp", - "stack": "cpp", - "xstring": "cpp", - "xutility": "cpp", - "xlocnum": "cpp", - "xtr1common": "cpp" - }, - "python.analysis.diagnosticSeverityOverrides": { - "reportInvalidTypeForm": "none" - } + "files.associations": { + "cctype": "cpp", + "clocale": "cpp", + "cmath": "cpp", + "cstdarg": "cpp", + "cstddef": "cpp", + "cstdio": "cpp", + "cstdlib": "cpp", + "cstring": "cpp", + "ctime": "cpp", + "cwchar": "cpp", + "cwctype": "cpp", + "array": "cpp", + "atomic": "cpp", + "strstream": "cpp", + "*.tcc": "cpp", + "bitset": "cpp", + "chrono": "cpp", + "cinttypes": "cpp", + "codecvt": "cpp", + "complex": "cpp", + "condition_variable": "cpp", + "cstdint": "cpp", + "deque": "cpp", + "forward_list": "cpp", + "list": "cpp", + "unordered_map": "cpp", + "unordered_set": "cpp", + "vector": "cpp", + "exception": "cpp", + "algorithm": "cpp", + "filesystem": "cpp", + "functional": "cpp", + "iterator": "cpp", + "map": "cpp", + "memory": "cpp", + "memory_resource": "cpp", + "numeric": "cpp", + "optional": "cpp", + "random": "cpp", + "ratio": "cpp", + "set": "cpp", + "string": "cpp", + "string_view": "cpp", + "system_error": "cpp", + "tuple": "cpp", + "type_traits": "cpp", + "utility": "cpp", + "fstream": "cpp", + "initializer_list": "cpp", + "iomanip": "cpp", + "iosfwd": "cpp", + "iostream": "cpp", + "istream": "cpp", + "limits": "cpp", + "mutex": "cpp", + "new": "cpp", + "ostream": "cpp", + "sstream": "cpp", + "stdexcept": "cpp", + "streambuf": "cpp", + "thread": "cpp", + "cfenv": "cpp", + "typeindex": "cpp", + "typeinfo": "cpp", + "valarray": "cpp", + "variant": "cpp", + "bit": "cpp", + "stack": "cpp", + "xstring": "cpp", + "xutility": "cpp", + "xlocnum": "cpp", + "xtr1common": "cpp" + }, + "python.analysis.diagnosticSeverityOverrides": { + "reportInvalidTypeForm": "none" + }, + "C_Cpp_Runner.cCompilerPath": "gcc", + "C_Cpp_Runner.cppCompilerPath": "g++", + "C_Cpp_Runner.debuggerPath": "gdb", + "C_Cpp_Runner.cStandard": "", + "C_Cpp_Runner.cppStandard": "", + "C_Cpp_Runner.msvcBatchPath": "C:/Program Files/Microsoft Visual Studio/VR_NR/Community/VC/Auxiliary/Build/vcvarsall.bat", + "C_Cpp_Runner.useMsvc": false, + "C_Cpp_Runner.warnings": [ + "-Wall", + "-Wextra", + "-Wpedantic", + "-Wshadow", + "-Wformat=2", + "-Wcast-align", + "-Wconversion", + "-Wsign-conversion", + "-Wnull-dereference" + ], + "C_Cpp_Runner.msvcWarnings": [ + "/W4", + "/permissive-", + "/w14242", + "/w14287", + "/w14296", + "/w14311", + "/w14826", + "/w44062", + "/w44242", + "/w14905", + "/w14906", + "/w14263", + "/w44265", + "/w14928" + ], + "C_Cpp_Runner.enableWarnings": true, + "C_Cpp_Runner.warningsAsError": false, + "C_Cpp_Runner.compilerArgs": [], + "C_Cpp_Runner.linkerArgs": [], + "C_Cpp_Runner.includePaths": [], + "C_Cpp_Runner.includeSearch": [ + "*", + "**/*" + ], + "C_Cpp_Runner.excludeSearch": [ + "**/build", + "**/build/**", + "**/.*", + "**/.*/**", + "**/.vscode", + "**/.vscode/**" + ], + "C_Cpp_Runner.useAddressSanitizer": false, + "C_Cpp_Runner.useUndefinedSanitizer": false, + "C_Cpp_Runner.useLeakSanitizer": false, + "C_Cpp_Runner.showCompilationTime": false, + "C_Cpp_Runner.useLinkTimeOptimization": false, + "C_Cpp_Runner.msvcSecureNoWarnings": false } \ No newline at end of file diff --git a/bxdf/.gitignore b/bxdf/.gitignore new file mode 100644 index 0000000..e7f4644 --- /dev/null +++ b/bxdf/.gitignore @@ -0,0 +1,3 @@ +*egg-info/ +dist/ +build/ \ No newline at end of file diff --git a/bxdf/medium.py b/bxdf/medium.py index 6fee29f..d5d76ef 100644 --- a/bxdf/medium.py +++ b/bxdf/medium.py @@ -23,6 +23,11 @@ class Medium_np: __type_mapping = {"hg": 0, "multi-hg": 1, "rayleigh": 2, "mie": 3, "transparent": -1} + + @staticmethod + def is_supported_type(_type: str): + return Medium_np.__type_mapping.get(_type, None) + def __init__(self, elem: xet.Element, is_world = False): """ Without spectral information, Rayleigh scattering here might not be physically-based diff --git a/bxdf/setup.py b/bxdf/setup.py new file mode 100644 index 0000000..8b222fe --- /dev/null +++ b/bxdf/setup.py @@ -0,0 +1,28 @@ +from pybind11.setup_helpers import Pybind11Extension +from setuptools import setup +import platform + +__version__ = "0.2.0" +cxx_std=11 + +ext_modules = [ + Pybind11Extension("vol_loader", + ["vol_loader/vol2numpy.cpp"], + # Example: passing in the version to the compiled code + define_macros = [('VERSION_INFO', __version__)], + extra_compile_args= ['-g', '-O3' if platform.system() == "Linux" else '-O2', '-fopenmp'], + ), +] + +setup( + name="vol_loader", + version=__version__, + author="Qianyue He", + description="Volume grid loader", + include_dirs="E:\\eigen-3.4.0", + long_description="", + ext_modules=ext_modules, + extras_require={"test": "pytest"}, + zip_safe=False, + python_requires=">=3.7", +) diff --git a/parsers/vol_loader/vol2numpy.cpp b/bxdf/vol_loader/vol2numpy.cpp similarity index 66% rename from parsers/vol_loader/vol2numpy.cpp rename to bxdf/vol_loader/vol2numpy.cpp index 41f4631..7a6747a 100644 --- a/parsers/vol_loader/vol2numpy.cpp +++ b/bxdf/vol_loader/vol2numpy.cpp @@ -6,6 +6,7 @@ * @date 2024-06-21 * @copyright Copyright (c) 2024 */ +#include #include #include #include @@ -20,8 +21,9 @@ struct VolumeData { int channels; std::vector data; - inline size_t size() const noexcept { - return xres * yres * zres; + inline size_t size(bool force_mono_color = false) const noexcept { + size_t ch = force_mono_color ? 1 : channels; + return ch * xres * yres * zres; } inline auto shape() const noexcept { @@ -71,46 +73,59 @@ bool readVolumeData(const std::string& filename, VolumeData& volume) { } // mitsuba3 vol to numpy -auto loadVol2Numpy(const std::string& filename) { +auto loadVol2Numpy(const std::string& filename, bool force_mono_color) { VolumeData volume; readVolumeData(filename, volume); - py::array_t vol_numpy(volume.size()); + py::array_t vol_numpy(volume.size(force_mono_color)); + float* const data_ptr = vol_numpy.mutable_data(0); if (volume.channels == 1) { + #pragma omp parallel for num_threads(4) for (int z = 0; z < volume.zres; ++z) { int zy_base = z * volume.yres; for (int y = 0; y < volume.yres; ++y) { int zyx_base = (zy_base + y) * volume.xres; for (int x = 0; x < volume.xres; ++x) { - vol_numpy[zyx_base + x] = volume.data[zyx_base + x]; + data_ptr[zyx_base + x] = volume.data[zyx_base + x]; } } } } else if (volume.channels == 3) { + #pragma omp parallel for num_threads(4) for (int z = 0; z < volume.zres; ++z) { int zy_base = z * volume.yres; for (int y = 0; y < volume.yres; ++y) { int zyx_base = (zy_base + y) * volume.xres; - for (int x = 0; x < volume.xres; ++x) { - int index = (zyx_base + x) * 3; - vol_numpy[index] = volume.data[index]; - vol_numpy[index + 1] = volume.data[index + 1]; - vol_numpy[index + 2] = volume.data[index + 2]; + if (!force_mono_color) { + for (int x = 0; x < volume.xres; ++x) { + int index = (zyx_base + x) * 3; + data_ptr[index] = volume.data[index]; + data_ptr[index + 1] = volume.data[index + 1]; + data_ptr[index + 2] = volume.data[index + 2]; + } + } else { + for (int x = 0; x < volume.xres; ++x) { + int index = zyx_base + x; + data_ptr[index] = volume.data[index * 3 + 1]; + } } } } } else { std::cerr << "Grid channel: <" << volume.channels << "> is not supported, supported channels: [1, 3]" << std::endl; } + + if (force_mono_color) + volume.channels = 1; - return std::tuple(vol_numpy, volume.shape()); + return std::tuple, std::tuple>(vol_numpy, volume.shape()); } PYBIND11_MODULE(vol_loader, m) { - m.doc() = "Volume grid (.vol / .vdb) loader (to numpy)\n"; m.def("vol_file_to_numpy", &loadVol2Numpy, "Load volume grid from mitsuba3 .vol file (return numpy)\n" - "Input: filename, input path of the file\n" + "Input: filename, [string] input path of the file\n" + "Input: force_mono_color, [bool] whether to extract only one channel from the volume data, default = False (True only for testing)\n" ); } \ No newline at end of file diff --git a/bxdf/volume.py b/bxdf/volume.py new file mode 100644 index 0000000..e637cfa --- /dev/null +++ b/bxdf/volume.py @@ -0,0 +1,464 @@ +""" + Grid Volume Scattering Medium (RGB / Float) + @author: Qianyue He + @date: 2024-6-22 +""" + +import os +import numpy as np +import taichi as ti +import taichi.math as tm +import xml.etree.ElementTree as xet + +from typing import Tuple +from taichi.math import vec2, vec3, vec4, mat3 +from bxdf.phase import PhaseFunction +from bxdf.medium import Medium_np +from la.cam_transform import delocalize_rotate +from parsers.general_parser import get, rgb_parse, transform_parse +from renderer.constants import ZERO_V3, ONES_V3 +from rich.console import Console +CONSOLE = Console(width = 128) + +vec3i = ti.types.vector(3, int) + +try: + from vol_loader import vol_file_to_numpy +except ImportError: + CONSOLE.log("[red]:skeleton: Error [/red]: vol_loader is not found, possible reason:") + CONSOLE.log(f"module not installed. Use 'python ./setup.py install --user' in {'./bxdf/'}") + raise ImportError("vol_loader not found") + +__all__ = ["GridVolume_np", "GridVolume"] + +FORCE_MONO_COLOR = False + +class GridVolume_np: + NONE = 0 + MONO = 1 + RGB = 2 + __type_mapping = {"none": NONE, "mono": MONO, "rgb": RGB} + def __init__(self, elem: xet.Element): + self.albedo = np.ones(3, np.float32) + self.phase_type_id = -1 + self.phase_type = "hg" + self.type_name = "none" + + self.type_id = GridVolume_np.NONE + self.xres = 0 + self.yres = 0 + self.zres = 0 + self.channel = 0 + + self.density_grid = None + + self.toWorld = None + self.scale = None + self.rotation = np.eye(3, dtype = np.float32) + self.offset = np.zeros(3, dtype = np.float32) + self.forward_t = None + self.density_scaling = np.ones(3, dtype = np.float32) + + # phase function + self.par = np.zeros(3, np.float32) + self.pdf = np.float32([1., 0., 0.]) + self.mono2rgb = False + + elem_to_query = { + "rgb": rgb_parse, + "float": lambda el: get(el, "value"), + "string": lambda el: self.setup_volume(get(el, "path", str)), + "transform": lambda el: self.assign_transform(*transform_parse(el)), + "bool": lambda el: el.get("value") in {"True", "true"} + } + if elem is not None: + type_name = elem.get("type") + if type_name in GridVolume_np.__type_mapping: + self.type_id = GridVolume_np.__type_mapping[type_name] + else: + CONSOLE.log(f"[red]Error :skull:[/red]: Volume '{elem.get('name')}' has unsupported type '{type_name}'") + raise NotImplementedError(f"GridVolume type '{type_name}' is not supported.") + self.type_name = type_name + + phase_type = elem.get("phase_type") + _phase_type_id = Medium_np.is_supported_type(phase_type) + if _phase_type_id is not None: + self.phase_type_id = _phase_type_id + else: + raise NotImplementedError(f"Phase function type '{phase_type}' is not supported.") + + for tag, query_func in elem_to_query.items(): + tag_elems = elem.findall(tag) + for tag_elem in tag_elems: + name = tag_elem.get("name") + if hasattr(self, name): + self.__setattr__(name, query_func(tag_elem)) + + if self.mono2rgb == True: + if self.channel == 3: + CONSOLE.log("Setting 'mono2rgb = True' is meanless for RGB volume. This setting is only meaningful when the volume is mono-chromatic.") + else: + CONSOLE.log("Mono-chromatic volume is converted to RGB volume.") + self.channel = 3 + self.type_id = GridVolume_np.RGB + self.density_grid = GridVolume_np.make_colorful_volume(self.density_grid, self.xres, self.yres, self.zres) + else: + self.density_grid = np.concatenate([self.density_grid, self.density_grid, self.density_grid], axis = -1) + + if self.scale is None: + self.scale = np.eye(3, dtype = np.float32) + else: + self.scale = np.diag(self.scale) + self.forward_t = self.rotation @ self.scale + if self.type_id == GridVolume_np.MONO: + self.density_grid *= self.density_scaling[0] + CONSOLE.log(f"Volume density scaled by {self.density_scaling[0]}.") + else: + self.density_grid *= self.density_scaling + CONSOLE.log(f"Volume density scaled by {self.density_scaling}.") + + def assign_transform(self, trans_r, trans_t, trans_s): + self.rotation = trans_r + self.offset = trans_t + self.scale = trans_s + + def setup_volume(self, path:str): + if not os.path.exists(path): + CONSOLE.log(f"[red]Error :skull:[/red]: {path} contains no valid volume file.") + raise RuntimeError("Volume file not found.") + density_grid, (self.xres, self.yres, self.zres, self.channel) = vol_file_to_numpy(path, FORCE_MONO_COLOR) + if FORCE_MONO_COLOR: + CONSOLE.log(f"[yellow]Warning[/yellow]: FORCE_MONO_COLOR is True. This only makes sense when we are testing the code.") + + density_grid = density_grid.reshape((self.zres, self.yres, self.xres, self.channel)) + return density_grid + + @staticmethod + def make_colorful_volume(density_grid: np.ndarray, xres: int, yres: int, zres: int): + density_grid = np.concatenate([ + density_grid, + density_grid, + density_grid + ], axis = -1) + + half_x = zres // 3 + grad_l = np.linspace(1, 0, half_x, dtype = np.float32) ** 0.65 + grad_r = np.linspace(0, 1, zres - half_x, dtype = np.float32) ** 0.6 + + left_half = np.zeros((half_x, 3), dtype = np.float32) + left_half[:, 0] = 1 - grad_l + left_half[:, 1] = 1 + left_half[:, 2] = 1 + + right_half = np.zeros((zres - half_x, 3), dtype = np.float32) + right_half[:, 0] = 1 + right_half[:, 1] = 1 + right_half[:, 2] = 1 - grad_r + + color_gradient = np.vstack((left_half, right_half)) + + return density_grid * color_gradient[:, None, None, :] + + def get_shape(self) -> Tuple[int, int, int]: + return (self.zres, self.yres, self.xres) + + def get_majorant(self, guard = 0.2, scale_ratio = 1.05): + maj = self.density_grid.max(axis = (0, 1, 2)) + maj_guard = np.mean(maj) * guard + maj = np.maximum(maj, maj_guard) + maj *= scale_ratio + return vec3(maj) if self.type_id == GridVolume_np.RGB else vec3([maj, maj, maj]) + + def export(self): + if self.type_id == GridVolume_np.NONE: + return GridVolume(_type = 0) + aabb_mini, aabb_maxi = self.get_aabb() + majorant = self.get_majorant() + CONSOLE.log(f"Majorant: {majorant}. PDF: {majorant / majorant.sum()}") + # the shape of density grid: (zres, yres, xres, channels) + return GridVolume( + _type = self.type_id, + albedo = vec3(self.albedo), + inv_T = mat3(np.linalg.inv(self.forward_t)), + trans = vec3(self.offset), + mini = vec3(aabb_mini), + maxi = vec3(aabb_maxi), + max_idxs = vec3i([self.xres - 1, self.yres - 1, self.zres - 1]), + majorant = majorant, + pdf = majorant / majorant.sum(), + ph = PhaseFunction(_type = self.phase_type_id, par = vec3(self.par), pdf = vec3(self.pdf)) + ) + + def local_to_world(self, point: np.ndarray) -> np.ndarray: + """ Take a point (shape (3) and transform it to the world space) """ + return point @ self.forward_t.T + self.offset + + def get_aabb(self): + x, y, z = self.xres, self.yres, self.zres + all_points = np.float32([ + [0, 0, 0], + [x, 0, 0], + [0, y, 0], + [x, y, 0], + [0, 0, z], + [x, 0, z], + [0, y, z], + [x, y, z] + ]) + world_point = self.local_to_world(all_points) + # conservative AABB + return world_point.min(axis = 0) - 0.01, world_point.max(axis = 0) + 0.01 + + + def __repr__(self): + aabb_min, aabb_max = self.get_aabb() + return f"" + + +@ti.func +def rgb_select(rgb: vec3, channel: int) -> float: + """ Get the value of a specified channel + This implement looks dumb, but I think it has its purpose: + Dynamic indexing will cause the local array get moved to global memory + i.e. the access speed will drop significantly + + I am not sure whether Taichi Lang has the problem or not, yet I know CUDA + has this problem + """ + result = 0.0 + if channel == 0: + result = rgb[0] + elif channel == 1: + result = rgb[1] + else: + result = rgb[2] + return result + +@ti.dataclass +class GridVolume: + """ Grid Volume Taichi End definition""" + _type: int + """ grid volume type: 0 (none), 1 (mono-chromatic), 2 (RGB) """ + albedo: vec3 + """ scattering albedo """ + inv_T: mat3 + """ Inverse transform matrix: $(R_rotation @ R_scaling)^{-1}$ """ + trans: vec3 + """ Translation vector (world frame) """ + mini: vec3 + """ Min bound (AABB) minimum coordinate """ + maxi: vec3 + """ Max bound (AABB) maximum coordinate """ + max_idxs: vec3i # max_indices + """ Maximum index for clamping, in case there's out-of-range access (xres - 1, yres - 1, zres - 1)""" + majorant: vec3 # max_values + """ (Scaled) Maximum value of sigma_t as majorant, for mono-chromatic, [0] is the value """ + pdf: vec3 # normalized majorant + """ PDF used for MIS (wavelength spectrum) """ + ph: PhaseFunction # phase function + """ Phase function for scattering sampling """ + + @ti.func + def is_scattering(self): # check whether the current medium is scattering medium + return self._type >= 1 + + @ti.func + def intersect_volume(self, ray_o: vec3, ray_d: vec3, max_t: float) -> vec2: + """ Get ray-volume intersection AABB near and far distance """ + near_far = vec2([0, max_t]) + inv_dir = 1.0 / ray_d + + t1s = (self.mini - ray_o) * inv_dir + t2s = (self.maxi - ray_o) * inv_dir + + tmin = ti.min(t1s, t2s) + tmax = ti.max(t1s, t2s) + + near_far[0] = ti.max(0, tmin.max()) + 1e-5 + near_far[1] = ti.min(max_t, tmax.min()) - 1e-5 + return near_far + + @ti.func + def transmittance(self, grid: ti.template(), ray_o: vec3, ray_d: vec3, thp: vec3, max_t: float) -> vec3: + transm = vec3([1, 1, 1]) + if self._type: + near_far = self.intersect_volume(ray_o, ray_d, max_t) + if near_far[0] < near_far[1] and near_far[1] > 0: + ray_o = self.inv_T @ (ray_o - self.trans) + ray_d = self.inv_T @ ray_d + if self._type: + transm = self.eval_tr_ratio_tracking_3d(grid, ray_o, ray_d, thp, near_far) + return transm + + @ti.func + def sample_mfp(self, grid: ti.template(), ray_o: vec3, ray_d: vec3, thp: vec3, max_t: float) -> vec4: + transm = vec4([1, 1, 1, -1]) + if self._type: + near_far = self.intersect_volume(ray_o, ray_d, max_t) + if near_far[0] < near_far[1] and near_far[1] > 0: + ray_o = self.inv_T @ (ray_o - self.trans) + ray_d = self.inv_T @ ray_d + if self._type: + transm = self.sample_distance_delta_tracking_3d(grid, ray_o, ray_d, thp, near_far) + return transm + + @ti.func + def density_lookup_3d(self, grid: ti.template(), index: vec3, u_offset: vec3) -> vec3: + """ Stochastic lookup of density (mono-chromatic volume) """ + idx = ti.cast(ti.floor(index + (u_offset - 0.5)), int) + val = ZERO_V3 + if (idx >= 0).all() and (idx <= self.max_idxs).all(): + val = grid[idx[2], idx[1], idx[0]] + return val + + @ti.func + def density_lookup_lerp_3d(self, grid: ti.template(), index: vec3, u_offset: vec3) -> vec3: + """ Stochastic lookup of density (mono-chromatic volume) """ + coord = index + (u_offset - 0.5) + idx = ti.cast(ti.floor(coord), int) + coord -= idx + val = ZERO_V3 + if (idx >= 0).all() and (idx <= self.max_idxs - 1).all(): + v1 = tm.mix(grid[idx[2], idx[1], idx[0]], grid[idx[2], idx[1], idx[0] + 1], coord[0]) + v2 = tm.mix(grid[idx[2] + 1, idx[1], idx[0]], grid[idx[2] + 1, idx[1], idx[0] + 1], coord[0]) + v3 = tm.mix(grid[idx[2], idx[1] + 1, idx[0]], grid[idx[2], idx[1] + 1, idx[0] + 1], coord[0]) + v4 = tm.mix(grid[idx[2] + 1, idx[1] + 1, idx[0]], grid[idx[2] + 1, idx[1] + 1, idx[0] + 1], coord[0]) + + v1 = tm.mix(v1, v2, coord[2]) + v3 = tm.mix(v3, v4, coord[2]) + + val = tm.mix(v1, v3, coord[1]) + return val + + @ti.func + def sample_new_rays(self, incid: vec3): + ret_spec = vec3([1, 1, 1]) + ret_dir = incid + ret_pdf = 1.0 + if self.is_scattering(): # medium interaction - evaluate phase function (currently output a float) + local_new_dir, ret_pdf = self.ph.sample_p(incid) # local frame ray_dir should be transformed + ret_dir, _ = delocalize_rotate(incid, local_new_dir) + ret_spec *= ret_pdf + return ret_dir, ret_spec, ret_pdf + + @ti.func + def sample_distance_delta_tracking_3d(self, + grid: ti.template(), ray_ol: vec3, ray_dl: vec3, thp: vec3, near_far: vec2 + ) -> vec4: + result = vec4([1, 1, 1, -1]) + if self._type: + # Step 1: choose one element according to the majorant + pdfs = thp * self.pdf + pdfs /= pdfs.sum() + Tr = 1.0 + pdf = 1.0 + albedo = 0.0 + maj = 0.0 + channel = 0 + val = ti.random(float) + + # avoid dynamic indexing on GPU (and array might be moved from local registers to global memory) + if val <= pdfs[0]: + albedo = self.albedo[0] + maj = self.majorant[0] + pdf = pdfs[0] + channel = 0 + elif val <= pdfs[0] + pdfs[1]: + albedo = self.albedo[1] + maj = self.majorant[1] + pdf = pdfs[1] + channel = 1 + else: + albedo = self.albedo[2] + maj = self.majorant[2] + pdf = pdfs[2] + channel = 2 + inv_maj = 1.0 / maj + + t = near_far[0] - ti.log(1.0 - ti.random(float)) * inv_maj + while t < near_far[1]: + d = self.density_lookup_3d(grid, ray_ol + t * ray_dl, vec3([ + ti.random(float), ti.random(float), ti.random(float) + ])) + # Scatter upon real collision + n_s = rgb_select(d, channel) + if ti.random(float) < n_s * inv_maj: + Tr *= albedo + result[3] = t + break + + t -= ti.log(1.0 - ti.random(float)) * inv_maj + if self._type == GridVolume_np.RGB: + if channel == 0: + result[:3] = vec3([Tr / pdf, 0, 0]) + elif channel == 1: + result[:3] = vec3([0, Tr / pdf, 0]) + else: + result[:3] = vec3([0, 0, Tr / pdf]) + else: + result[:3] = vec3([Tr, Tr, Tr]) + return result + + @ti.func + def eval_tr_ratio_tracking_3d(self, + grid: ti.template(), ray_ol: vec3, ray_dl: vec3, thp: vec3, near_far: vec2 + ) -> vec3: + transm = ONES_V3 + if self._type: + # Step 1: choose one element according to the majorant + pdfs = thp * self.pdf + pdfs /= pdfs.sum() + Tr = 1.0 + pdf = 1.0 + maj = 0.0 + channel = 0 + val = ti.random(float) + + # avoid dynamic indexing on GPU (and array might be moved from local registers to global memory) + if val <= pdfs[0]: + maj = self.majorant[0] + pdf = pdfs[0] + channel = 0 + elif val <= pdfs[0] + pdfs[1]: + maj = self.majorant[1] + pdf = pdfs[1] + channel = 1 + else: + maj = self.majorant[2] + pdf = pdfs[2] + channel = 2 + inv_maj = 1.0 / maj + + t = near_far[0] + while True: + # problem with coordinates + t -= ti.log(1.0 - ti.random(float)) * inv_maj + + if t >= near_far[1]: break + # for mono-chromatic medium, this is 1 + d = self.density_lookup_3d(grid, ray_ol + t * ray_dl, vec3([ + ti.random(float), ti.random(float), ti.random(float) + ])) + + n_s = rgb_select(d, channel) + Tr *= ti.max(0.0, 1.0 - n_s * inv_maj) + + # Russian Roulette + if Tr < 0.1: + if ti.random(float) >= Tr: + Tr = 0.0 + break + Tr = 1.0 + if self._type == GridVolume_np.RGB: + if channel == 0: + transm = vec3([Tr / pdf, 0, 0]) + elif channel == 1: + transm = vec3([0, Tr / pdf, 0]) + else: + transm = vec3([0, 0, Tr / pdf]) + else: + transm = vec3([Tr, Tr, Tr]) + return transm + \ No newline at end of file diff --git a/emitters/abtract_source.py b/emitters/abtract_source.py index 761dda0..720733c 100644 --- a/emitters/abtract_source.py +++ b/emitters/abtract_source.py @@ -107,8 +107,6 @@ def sample_hit( if ti.static(HEMISPHERE_SAMPLE_SPHERE): to_hit = (hit_pos - center).normalized() - # the pdf here can be viewed as being both both sa & area measure - # since for a unit sphere, different unit solid angle extends to the same amount of area local_dir, pdf = uniform_sphere() normal, _ = delocalize_rotate(to_hit, local_dir) ret_pos = center + normal * radius diff --git a/parsers/general_parser.py b/parsers/general_parser.py index 42bffc1..0df0f71 100644 --- a/parsers/general_parser.py +++ b/parsers/general_parser.py @@ -53,32 +53,34 @@ def vec3d_parse(elem: xet.Element): # for positions, implicit float -> vec3 conversion is not allowed return parse_str(elem.get("value"), no_else_branch = True) -def transform_parse(transform_elem: xet.Element) -> Tuple[Arr, Arr]: +def transform_parse(transform_elem: xet.Element) -> Tuple[Arr, Arr, Arr]: """ Note that: extrinsic rotation is not supported, meaning that we can only rotate around the object centroid, which is [intrinsic rotation]. Yet, extrinsic rotation can be composed by intrinsic rotation and translation """ - trans_r, trans_t = None, None + trans_r, trans_t, trans_s = None, None, None for child in transform_elem: if child.tag == "translate": trans_t = np.float32([get(child, "x"), get(child, "y"), get(child, "z")]) elif child.tag == "rotate": - rot_type = child.get("type") + rot_type = child.get("type", "euler") if rot_type == "euler": r_angle = get(child, "r") # roll p_angle = get(child, "p") # pitch y_angle = get(child, "y") # yaw trans_r = Rot.from_euler("zxy", (r_angle, p_angle, y_angle), degrees = True).as_matrix() elif rot_type == "quaternion": - trans_r = Rot.from_quat([get(child, "x"), get(child, "y"), get(child, "z"), get(child, "w")]) + trans_r = Rot.from_quat([get(child, "x"), get(child, "y"), get(child, "z"), get(child, "w")]).as_matrix() elif rot_type == "angle-axis": axis: Arr = np.float32([get(child, "x"), get(child, "y"), get(child, "z")]) axis /= np.linalg.norm(axis) * get(child, "angle") / 180. * np.pi - trans_r = Rot.from_rotvec(axis) + trans_r = Rot.from_rotvec(axis).as_matrix() else: raise ValueError(f"Unsupported rotation representation '{rot_type}'") + elif child.tag == "scale": + trans_s = np.float32([get(child, "x"), get(child, "y"), get(child, "z")]) elif child.tag.lower() == "lookat": target_point = parse_str(child.get("target")) origin_point = parse_str(child.get("origin")) @@ -93,7 +95,7 @@ def transform_parse(transform_elem: xet.Element) -> Tuple[Arr, Arr]: raise ValueError(f"Unsupported transformation representation '{child.tag}'") # Note that, trans_r (rotation) is defualt to be intrinsic (apply under the centroid coordinate) # Therefore, do not use trans_r unless you know how to correctly transform objects with it - return trans_r, trans_t # trans_t and trans_r could be None, if is not defined in the object + return trans_r, trans_t, trans_s # trans_t, trans_r and trans_s could be None, if is not defined in the object def parse_sphere_element(elem: xet.Element): sphere_info = np.zeros((1, 2, 3), np.float32) diff --git a/parsers/obj_loader.py b/parsers/obj_loader.py index 0776dc5..54e63bd 100644 --- a/parsers/obj_loader.py +++ b/parsers/obj_loader.py @@ -92,13 +92,24 @@ def calculate_surface_area(meshes: Arr, _type = 0): area_sum = 4. * np.pi * radius ** 2 return area_sum -def apply_transform(meshes: Arr, normals: Arr, trans_r: Arr, trans_t: Arr) -> Arr: +def is_uniform_scaling(scale: Arr) -> bool: + if scale[0] != scale[1] or scale[0] != scale[2]: + return False + return True + +def apply_transform(meshes: Arr, normals: Arr, trans_r: Arr, trans_t: Arr, trans_s: Arr) -> Arr: """ - input normals are of shape (N, 3) - input meshes are of shape (N, 3, 3), and for the last two dims - 3(front): entry index for the vertices of triangles - 3(back): entry index for (x, y, z) """ + if trans_s is not None: + if not is_uniform_scaling(trans_s): + CONSOLE.log("Warning: scaling for meshes should be uniform, otherwise normals should be re-computed.") + CONSOLE.log(f"Scaling factor is adjusted to be: {trans_s[0]}") + trans_s[1] = trans_s[0] + trans_s[2] = trans_s[0] if trans_r is not None: center = meshes.mean(axis = 1).mean(axis = 0) meshes -= center # decentralize diff --git a/parsers/xml_parser.py b/parsers/xml_parser.py index 288619a..3bb5d47 100644 --- a/parsers/xml_parser.py +++ b/parsers/xml_parser.py @@ -86,6 +86,10 @@ def parse_emitters(em_elem: list): raise ValueError(f"Source [{emitter_type}] is not implemented.") return sources, source_id_dict +def parse_volume(directory: str, vol_list: List[xet.Element]): + + pass + def parse_wavefront( directory: str, obj_list: List[xet.Element], bsdf_dict: dict, emitter_dict: dict, texture_dict: List[Texture_np]) -> List[Arr]: @@ -111,8 +115,8 @@ def parse_wavefront( meshes, normals, vns, uvs = extract_obj_info(os.path.join(directory, filepath_child.get("value"))) transform_child = elem.find("transform") if transform_child is not None: - trans_r, trans_t = transform_parse(transform_child) - meshes, normals = apply_transform(meshes, normals, trans_r, trans_t) + trans_r, trans_t, trans_s = transform_parse(transform_child) + meshes, normals = apply_transform(meshes, normals, trans_r, trans_t, trans_s) if vns is not None: has_vertex_normal = True else: # CURRENTLY, only sphere is supported @@ -255,6 +259,10 @@ def scene_parsing(directory: str, file: str): shape_nodes = root_node.findall("shape") sensor_node = root_node.find("sensor") world_node = root_node.find("world") + volume_node = root_node.findall("volume") + if len(volume_node) > 1: + CONSOLE.log(f"There are {len(volume_node)} in total, only the first volume will be kept.") + volume_node = volume_node[:1] assert(sensor_node) emitter_configs, \ emitter_dict = parse_emitters(emitter_nodes) @@ -268,7 +276,7 @@ def scene_parsing(directory: str, file: str): Each mapping might needs a different packing, therefore we need different image packaging and texture info block For normal mapping, non-bidirectional renderers will be simple but not for BDPT roughness is of lower priority - - [ ] Speed up python->taichi conversion + - [x] Speed up python->taichi conversion """ array_info, all_objs, area_lut, has_vertex_normal \ = parse_wavefront(directory, shape_nodes, bsdf_dict, emitter_dict, textures) @@ -276,6 +284,7 @@ def scene_parsing(directory: str, file: str): configs['world'] = parse_world(world_node) configs['packed_textures'] = teximgs configs['has_vertex_normal'] = has_vertex_normal + configs['volume'] = volume_node emitter_configs = update_emitter_config(emitter_configs, area_lut) return emitter_configs, array_info, all_objs, configs diff --git a/renderer/bdpt.py b/renderer/bdpt.py index 99f3931..93006b8 100644 --- a/renderer/bdpt.py +++ b/renderer/bdpt.py @@ -240,7 +240,7 @@ def random_walk( # Step 2: check for mean free path sampling # Calculate mfp, path_beta = transmittance / PDF - is_mi, it.min_depth, path_beta = self.sample_mfp(it.obj_id, in_free_space, it.min_depth) + is_mi, it.min_depth, path_beta = self.sample_mfp(ray_o, ray_d, throughput, it.obj_id, in_free_space, it.min_depth) if it.obj_id < 0 and not is_mi: break # exiting world bound throughput *= path_beta # attenuate first if throughput.max() < 5e-5: break @@ -392,7 +392,7 @@ def connect_path(self, i: int, j: int, sid: int, tid: int, decomp: int): le = cam_v.beta * fr_cam * fr_lit * lit_v.beta / (depth * depth) ret_time = lit_v.time + cam_v.time if le.max() > 0 and calc_transmittance == True: - tr, track_depth = self.track_ray(connect_dir, track_pos, depth) + tr, track_depth = self.track_ray(connect_dir, track_pos, ONES_V3, depth) le *= tr ret_time += track_depth weight = 0. diff --git a/renderer/constants.py b/renderer/constants.py index 684c513..9b9f711 100644 --- a/renderer/constants.py +++ b/renderer/constants.py @@ -34,6 +34,7 @@ ZERO_V3 = vec3([0, 0, 0]) ONES_V3 = vec3([1, 1, 1]) +EPSL_V3 = vec3([1e-4, 1e-4, 1e-4]) INVALID = vec3([-1, -1, -1]) AXIS_X = vec3([1, 0, 0]) AXIS_Y = vec3([0, 1, 0]) diff --git a/renderer/vanilla_renderer.py b/renderer/vanilla_renderer.py index b8a4606..007d854 100644 --- a/renderer/vanilla_renderer.py +++ b/renderer/vanilla_renderer.py @@ -45,13 +45,14 @@ def render(self, _t_start: int, _t_end: int, _s_start: int, _s_end: int, _a: int color = vec3([0, 0, 0]) contribution = vec3([1, 1, 1]) emission_weight = 1.0 - for _i in range(self.max_bounce): + for bounce in range(self.max_bounce): if it.is_ray_not_hit(): break # nothing is hit, break if ti.static(self.use_rr): # Simple Russian Roullete ray termination max_value = contribution.max() - if ti.random(float) > max_value: break - else: contribution *= 1. / (max_value + 1e-7) # unbiased calculation + if max_value < self.rr_threshold and bounce >= self.rr_bounce_th: + if ti.random(float) > max_value: break + else: contribution *= 1. / (max_value + 1e-7) # unbiased calculation else: if contribution.max() < 1e-4: break # contribution too small, break hit_point = ray_d * it.min_depth + ray_o diff --git a/renderer/vpt.py b/renderer/vpt.py index 348347d..7527da5 100644 --- a/renderer/vpt.py +++ b/renderer/vpt.py @@ -13,6 +13,7 @@ from la.cam_transform import * from tracer.path_tracer import PathTracer from emitters.abtract_source import LightSource +from bxdf.volume import GridVolume_np from parsers.obj_desc import ObjDescriptor from sampler.general_sampling import balance_heuristic @@ -34,6 +35,20 @@ def __init__(self, objects: List[ObjDescriptor], prop: dict, bvh_delay: bool = False ): super().__init__(emitters, array_info, objects, prop, bvh_delay) + self.has_volume = False + if prop["volume"]: + CONSOLE.log("Loading Grid Volume ...") + volume = GridVolume_np(prop["volume"][0]) + self.density_grid = ti.Vector.field(3, float, volume.get_shape()) + if volume.type_id not in{ GridVolume_np.MONO, GridVolume_np.RGB}: + CONSOLE.log(f"Volume type '{volume.type_name}' is not supported.") + raise RuntimeError("Unsupported volume type.") + self.has_volume = True + self.density_grid.from_numpy(volume.density_grid) + self.volume = volume.export() + CONSOLE.log(f"Grid volume loaded: {volume}") + else: + self.density_grid = ti.field(float, (1, 1, 1)) self.world_scattering = self.world.medium._type >= 0 @ti.func @@ -57,11 +72,11 @@ def non_null_surface(self, idx: int): return non_null @ti.func - def sample_mfp(self, idx: int, in_free_space: int, depth: float): + def sample_mfp(self, ray_o: vec3, ray_d: vec3, thp: vec3, idx: int, in_free_space: int, depth: float): """ Mean free path sampling, returns: - whether medium is a valid scattering medium / mean free path """ - is_mi = False + is_mi = 0 mfp = depth beta = vec3([1., 1., 1.]) # whether the world is valid for scattering: inside the free space and world has scattering medium @@ -75,10 +90,16 @@ def sample_mfp(self, idx: int, in_free_space: int, depth: float): elif not in_free_space: is_mi, mfp, beta = self.bsdf_field[idx].medium.sample_mfp(depth) # use medium to sample / calculate transmittance + if ti.static(self.has_volume): # grid volume event might override the world medium event (not physically based, but simple to implement) + result = self.volume.sample_mfp(self.density_grid, ray_o, ray_d, thp, depth) + if result[3] > 0: # in case grid volume is nested with world medium + is_mi = 2 + mfp = result[3] + beta = result[:3] return is_mi, mfp, beta @ti.func - def track_ray(self, ray, start_p, depth): + def track_ray(self, cur_ray: vec3, cur_point: vec3, thp: vec3, depth: float): """ For medium interaction, check if the path to one point is not blocked (by non-null surface object) And also we need to calculate the attenuation along the path, e.g.: if the ray passes through @@ -86,9 +107,10 @@ def track_ray(self, ray, start_p, depth): FIXME: the speed of this method should be boosted """ tr = vec3([1., 1., 1.]) + if ti.static(self.has_volume): # grid volume transmittance + tr = self.volume.transmittance(self.density_grid, cur_point, cur_ray, thp, depth) + in_free_space = True - cur_point = start_p - cur_ray = ray acc_depth = 0.0 for _i in range(7): # maximum tracking depth = 7 (corresponding to at most 2 clouds of smoke) it = self.ray_intersect(cur_ray, cur_point, depth) @@ -133,19 +155,24 @@ def render(self, _t_start: int, _t_end: int, _s_start: int, _s_end: int, _a: int color = vec3([0, 0, 0]) throughput = vec3([1, 1, 1]) emission_weight = 1.0 - + in_free_space = True bounce = 0 while True: # for _i in range(self.max_bounce): - # Step 1: ray termination test - Only RR termination is allowed - max_value = throughput.max() - if ti.random(float) > max_value: break - else: throughput *= 1. / ti.max(max_value, 1e-7) # unbiased calculation + # Step 1: ray termination test + if ti.static(self.use_rr): + # Simple Russian Roullete ray termination + max_value = throughput.max() + if max_value < self.rr_threshold and bounce >= self.rr_bounce_th: + if ti.random(float) > max_value: break + else: throughput *= 1. / (max_value + 1e-7) # unbiased calculation + else: + if throughput.max() < 1e-5: break # contribution too small, break # Step 2: ray intersection it = self.ray_intersect(ray_d, ray_o) if it.obj_id < 0: - if not self.world_scattering: break # nothing is hit, break + if ti.static(not self.world_scattering and not self.has_volume): break # nothing is hit, break else: # the world is filled with scattering medium it.min_depth = self.world_bound_time(ray_o, ray_d) in_free_space = True @@ -154,8 +181,9 @@ def render(self, _t_start: int, _t_end: int, _s_start: int, _s_end: int, _a: int in_free_space = tm.dot(it.n_g, ray_d) < 0 # Step 3: check for mean free path sampling # Calculate mfp, path_beta = transmittance / PDF - is_mi, it.min_depth, path_beta = self.sample_mfp(it.obj_id, in_free_space, it.min_depth) + is_mi, it.min_depth, path_beta = self.sample_mfp(ray_o, ray_d, throughput, it.obj_id, in_free_space, it.min_depth) if it.obj_id < 0 and not is_mi: break # exiting world bound + hit_point = ray_d * it.min_depth + ray_o throughput *= path_beta # attenuate first if not is_mi and not self.non_null_surface(it.obj_id): @@ -180,7 +208,7 @@ def render(self, _t_start: int, _t_end: int, _s_start: int, _s_end: int, _a: int to_emitter = emit_pos - hit_point emitter_d = to_emitter.norm() light_dir = to_emitter / emitter_d - tr, _ = self.track_ray(light_dir, hit_point, emitter_d) + tr, _ = self.track_ray(light_dir, hit_point, throughput, emitter_d) shadow_int *= tr direct_spec = self.eval(it, ray_d, light_dir, is_mi, in_free_space) else: # the only situation for being invalid, is when there is only one source and the ray hit the source diff --git a/scenes/cbox/cbox-rgbvol.xml b/scenes/cbox/cbox-rgbvol.xml new file mode 100644 index 0000000..3f09345 --- /dev/null +++ b/scenes/cbox/cbox-rgbvol.xml @@ -0,0 +1,171 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/scenes/cbox/cbox-volgrid.xml b/scenes/cbox/cbox-volgrid.xml new file mode 100755 index 0000000..e0e42f9 --- /dev/null +++ b/scenes/cbox/cbox-volgrid.xml @@ -0,0 +1,185 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/scenes/csphere/balls-mono.xml b/scenes/csphere/balls-mono.xml index 8a8a410..cdbd0f6 100644 --- a/scenes/csphere/balls-mono.xml +++ b/scenes/csphere/balls-mono.xml @@ -78,7 +78,7 @@ - + diff --git a/scenes/meshes/cornell/cbox_blue_light.obj b/scenes/meshes/cornell/cbox_blue_light.obj new file mode 100644 index 0000000..276088e --- /dev/null +++ b/scenes/meshes/cornell/cbox_blue_light.obj @@ -0,0 +1,15 @@ +# Blender v3.0.1 OBJ File: '' +# www.blender.org +o Plane +v 0.042202 1.211291 5.535800 +v 5.470666 1.211291 5.535800 +v 0.042202 0.884762 5.535800 +v 5.470666 0.884762 5.535800 +vt 1.000000 0.000000 +vt 0.000000 1.000000 +vt 0.000000 0.000000 +vt 1.000000 1.000000 +vn 0.0000 0.0000 -1.0000 +s off +f 2/1/1 3/2/1 1/3/1 +f 2/1/1 4/4/1 3/2/1 diff --git a/scenes/meshes/cornell/cbox_green_light.obj b/scenes/meshes/cornell/cbox_green_light.obj new file mode 100644 index 0000000..184111f --- /dev/null +++ b/scenes/meshes/cornell/cbox_green_light.obj @@ -0,0 +1,15 @@ +# Blender v3.0.1 OBJ File: '' +# www.blender.org +o Plane.001 +v 0.042202 2.768821 5.535800 +v 5.470666 2.768821 5.535800 +v 0.042202 2.442292 5.535800 +v 5.470666 2.442292 5.535800 +vt 1.000000 0.000000 +vt 0.000000 1.000000 +vt 0.000000 0.000000 +vt 1.000000 1.000000 +vn 0.0000 0.0000 -1.0000 +s off +f 2/1/1 3/2/1 1/3/1 +f 2/1/1 4/4/1 3/2/1 diff --git a/scenes/meshes/cornell/cbox_red_light.obj b/scenes/meshes/cornell/cbox_red_light.obj new file mode 100644 index 0000000..67c8223 --- /dev/null +++ b/scenes/meshes/cornell/cbox_red_light.obj @@ -0,0 +1,15 @@ +# Blender v3.0.1 OBJ File: '' +# www.blender.org +o Plane.002 +v 0.042202 4.489385 5.535800 +v 5.470666 4.489385 5.535800 +v 0.042202 4.162856 5.535800 +v 5.470666 4.162856 5.535800 +vt 1.000000 0.000000 +vt 0.000000 1.000000 +vt 0.000000 0.000000 +vt 1.000000 1.000000 +vn 0.0000 0.0000 -1.0000 +s off +f 2/1/1 3/2/1 1/3/1 +f 2/1/1 4/4/1 3/2/1 diff --git a/tracer/path_tracer.py b/tracer/path_tracer.py index da2168f..5cb70c3 100644 --- a/tracer/path_tracer.py +++ b/tracer/path_tracer.py @@ -9,7 +9,6 @@ import sys sys.path.append("..") -import tqdm import numpy as np import taichi as ti import taichi.math as tm @@ -22,6 +21,7 @@ from bxdf.brdf import BRDF from bxdf.bsdf import BSDF, BSDF_np +from bxdf.volume import GridVolume_np from bxdf.texture import Texture, Texture_np from parsers.opts import get_options from parsers.obj_desc import ObjDescriptor @@ -31,7 +31,7 @@ from sampler.general_sampling import * from utils.tools import TicToc from tracer.interaction import Interaction -from tracer.ti_bvh import LinearBVH, LinearNode, export_python_bvh +from tracer.ti_bvh import LinearBVH, LinearNode from rich.console import Console CONSOLE = Console(width = 128) @@ -64,6 +64,8 @@ def __init__(self, self.stratified_sample = prop['stratified_sampling'] # whether to use stratified sampling self.use_mis = prop['use_mis'] # whether to use multiple importance sampling self.num_shadow_ray = prop['num_shadow_ray'] # number of shadow samples to trace + self.rr_threshold = prop.get('rr_threshold', 0.1) # threshold of employing RR + self.rr_bounce_th = prop.get('rr_bounce_th', 4) # minimum number of bounces to start RR self.brdf_two_sides = prop.get('brdf_two_sides', False) if self.num_shadow_ray > 0: @@ -90,6 +92,8 @@ def __init__(self, # maybe we should opt for environment map (env lighting) self.roughness_map = Texture.field() # TODO: this is useless for now + self.volume = GridVolume_np(None).export() + ti.root.dense(ti.i, self.src_num).place(self.src_field) # Light source Taichi storage self.obj_nodes = ti.root.bitmasked(ti.i, self.num_objects) self.obj_nodes.place(self.brdf_field) # BRDF Taichi storage @@ -433,10 +437,13 @@ def sample_new_ray(self, ret_pdf = 1.0 is_specular = False if is_mi: - if in_free_space: # sample world medium - ret_dir, ret_spec, ret_pdf = self.world.medium.sample_new_rays(incid) - else: # sample object medium - ret_dir, ret_spec, ret_pdf = self.bsdf_field[it.obj_id].medium.sample_new_rays(incid) + if is_mi > 1: # well this is..., I know it is ugly but... + ret_dir, ret_spec, ret_pdf = self.volume.sample_new_rays(incid) + else: + if in_free_space: # sample world medium + ret_dir, ret_spec, ret_pdf = self.world.medium.sample_new_rays(incid) + else: # sample object medium + ret_dir, ret_spec, ret_pdf = self.bsdf_field[it.obj_id].medium.sample_new_rays(incid) else: # surface sampling if ti.is_active(self.obj_nodes, it.obj_id): # active means the object is attached to BRDF if ti.static(self.brdf_two_sides): diff --git a/tracer/setup.py b/tracer/setup.py index ecdae95..0646bf6 100644 --- a/tracer/setup.py +++ b/tracer/setup.py @@ -1,5 +1,6 @@ from pybind11.setup_helpers import Pybind11Extension from setuptools import setup +import platform __version__ = "0.1.0" cxx_std=11 @@ -9,7 +10,7 @@ ["bvh/bvh.cpp"], # Example: passing in the version to the compiled code define_macros = [('VERSION_INFO', __version__)], - extra_compile_args= ['-g', '-O3'], + extra_compile_args= ['-g', '-O3' if platform.system() == "Linux" else '-O2'], ), ] @@ -18,7 +19,7 @@ version=__version__, author="Qianyue He", description="BVH constructed via C++ backend", - include_dirs="/usr/include/eigen3/", + include_dirs="E:\\eigen-3.4.0", long_description="", ext_modules=ext_modules, extras_require={"test": "pytest"},