Skip to content

Adding SDF format to kernel. #324

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
meet-brad-ch wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Adding SDF format to kernel. #324

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
88 changes: 88 additions & 0 deletions src/kernel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2383,6 +2383,94 @@ string opencl_c_container() { return R( // ########################## begin of O
}
} // voxelize_mesh()

)+R(kernel void voxelize_sdf)+"("+R(global uchar* flags, const uchar flag, const global float* sdf_data, const global float* sdf_params)+") {"+R( // voxelize_sdf()
const uxx n = get_global_id(0);
if(n>=(uxx)def_N) return;

// Unpack parameters
const uint sdf_Nx = as_uint(sdf_params[0]);
const uint sdf_Ny = as_uint(sdf_params[1]);
const uint sdf_Nz = as_uint(sdf_params[2]);
const float sdf_cell_size = sdf_params[3];
const float3 sdf_origin = (float3)(sdf_params[4], sdf_params[5], sdf_params[6]);
const float3 center_lbm = (float3)(sdf_params[7], sdf_params[8], sdf_params[9]);
const float world_to_lbm = sdf_params[10];
// Inverse rotation matrix (row-major)
const float3 rot_row0 = (float3)(sdf_params[11], sdf_params[12], sdf_params[13]);
const float3 rot_row1 = (float3)(sdf_params[14], sdf_params[15], sdf_params[16]);
const float3 rot_row2 = (float3)(sdf_params[17], sdf_params[18], sdf_params[19]);

// Get LBM coordinates
const uint3 xyz = coordinates(n);
const float3 lbm_pos = (float3)(xyz.x + 0.5f, xyz.y + 0.5f, xyz.z + 0.5f);

// Transform LBM -> World -> SDF coordinates
// 1. Compute offset from LBM center
const float3 offset_lbm = (lbm_pos - center_lbm) / world_to_lbm;
// 2. Apply inverse rotation to get offset in SDF world space
const float3 offset_world = (float3)(dot(rot_row0, offset_lbm), dot(rot_row1, offset_lbm), dot(rot_row2, offset_lbm));
// 3. Add SDF world center to get world position
const float3 sdf_world_size = (float3)(sdf_Nx, sdf_Ny, sdf_Nz) * sdf_cell_size;
const float3 sdf_world_center = sdf_origin + 0.5f * sdf_world_size;
const float3 world_pos = offset_world + sdf_world_center;
// 4. Convert to SDF grid coordinates
const float3 sdf_coord = (world_pos - sdf_origin) / sdf_cell_size;

// Check bounds
if(sdf_coord.x < 0.0f || sdf_coord.x >= (float)(sdf_Nx-1u) ||
sdf_coord.y < 0.0f || sdf_coord.y >= (float)(sdf_Ny-1u) ||
sdf_coord.z < 0.0f || sdf_coord.z >= (float)(sdf_Nz-1u)) {
return; // Outside SDF bounds
}

// Trilinear interpolation
const uint x0 = (uint)sdf_coord.x;
const uint y0 = (uint)sdf_coord.y;
const uint z0 = (uint)sdf_coord.z;
const uint x1 = min(x0 + 1u, sdf_Nx - 1u);
const uint y1 = min(y0 + 1u, sdf_Ny - 1u);
const uint z1 = min(z0 + 1u, sdf_Nz - 1u);

const float fx = sdf_coord.x - (float)x0;
const float fy = sdf_coord.y - (float)y0;
const float fz = sdf_coord.z - (float)z0;

// Sample SDF at 8 corners (SDF uses x + Nx*(y + Ny*z) indexing)
const ulong sdf_index_000 = (ulong)x0 + (ulong)sdf_Nx * ((ulong)y0 + (ulong)sdf_Ny * (ulong)z0);
const ulong sdf_index_001 = (ulong)x0 + (ulong)sdf_Nx * ((ulong)y0 + (ulong)sdf_Ny * (ulong)z1);
const ulong sdf_index_010 = (ulong)x0 + (ulong)sdf_Nx * ((ulong)y1 + (ulong)sdf_Ny * (ulong)z0);
const ulong sdf_index_011 = (ulong)x0 + (ulong)sdf_Nx * ((ulong)y1 + (ulong)sdf_Ny * (ulong)z1);
const ulong sdf_index_100 = (ulong)x1 + (ulong)sdf_Nx * ((ulong)y0 + (ulong)sdf_Ny * (ulong)z0);
const ulong sdf_index_101 = (ulong)x1 + (ulong)sdf_Nx * ((ulong)y0 + (ulong)sdf_Ny * (ulong)z1);
const ulong sdf_index_110 = (ulong)x1 + (ulong)sdf_Nx * ((ulong)y1 + (ulong)sdf_Ny * (ulong)z0);
const ulong sdf_index_111 = (ulong)x1 + (ulong)sdf_Nx * ((ulong)y1 + (ulong)sdf_Ny * (ulong)z1);

const float c000 = sdf_data[sdf_index_000];
const float c001 = sdf_data[sdf_index_001];
const float c010 = sdf_data[sdf_index_010];
const float c011 = sdf_data[sdf_index_011];
const float c100 = sdf_data[sdf_index_100];
const float c101 = sdf_data[sdf_index_101];
const float c110 = sdf_data[sdf_index_110];
const float c111 = sdf_data[sdf_index_111];

// Trilinear interpolation
const float c00 = c000 * (1.0f - fx) + c100 * fx;
const float c01 = c001 * (1.0f - fx) + c101 * fx;
const float c10 = c010 * (1.0f - fx) + c110 * fx;
const float c11 = c011 * (1.0f - fx) + c111 * fx;

const float c0 = c00 * (1.0f - fy) + c10 * fy;
const float c1 = c01 * (1.0f - fy) + c11 * fy;

const float sdf_value = c0 * (1.0f - fz) + c1 * fz;

// Mark as solid if inside (negative SDF)
if(sdf_value <= 0.0f) {
flags[n] = (flags[n] & ~TYPE_BO) | flag;
}
} // voxelize_sdf()

)+R(kernel void unvoxelize_mesh(global uchar* flags, const uchar flag, float x0, float y0, float z0, float x1, float y1, float z1) { // remove voxelized triangle mesh
const uxx n = get_global_id(0);
const float3 p = position(coordinates(n))+(float3)(0.5f*(float)((int)def_Nx+2*def_Ox)-0.5f, 0.5f*(float)((int)def_Ny+2*def_Oy)-0.5f, 0.5f*(float)((int)def_Nz+2*def_Oz)-0.5f);
Expand Down
73 changes: 73 additions & 0 deletions src/lbm.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -325,6 +325,43 @@ void LBM_Domain::voxelize_mesh_on_device(const Mesh* mesh, const uchar flag, con
bounding_box_and_velocity.write_to_device();
kernel_voxelize_mesh.run();
}
void LBM_Domain::voxelize_sdf_on_device(const SDF* sdf, const float3& center, const float3x3& rotation, const float scale, const uchar flag) { // voxelize SDF on GPU
// Upload SDF data to GPU
const ulong sdf_size = (ulong)sdf->Nx * (ulong)sdf->Ny * (ulong)sdf->Nz;
Memory<float> sdf_data(device, sdf_size, 1u, sdf->data);

// Compute inverse rotation matrix (to transform LBM coords back to SDF space)
// For orthogonal matrices, inverse = transpose
const float3x3 rotation_inv = transpose(rotation);

// Pack parameters into array (similar to voxelize_mesh approach)
Memory<float> sdf_params(device, 20u);
const float3 sdf_world_size = sdf->get_world_size();
const float sdf_max_world_dim = fmax(fmax(sdf_world_size.x, sdf_world_size.y), sdf_world_size.z);
const float world_to_lbm = scale / sdf_max_world_dim;

sdf_params[0] = as_float(sdf->Nx); // pack uint as float
sdf_params[1] = as_float(sdf->Ny);
sdf_params[2] = as_float(sdf->Nz);
sdf_params[3] = sdf->cell_size;
sdf_params[4] = sdf->origin.x; // SDF origin
sdf_params[5] = sdf->origin.y;
sdf_params[6] = sdf->origin.z;
sdf_params[7] = center.x; // LBM center
sdf_params[8] = center.y;
sdf_params[9] = center.z;
sdf_params[10] = world_to_lbm;
// Inverse rotation matrix (row-major: row0, row1, row2)
sdf_params[11] = rotation_inv.xx; sdf_params[12] = rotation_inv.xy; sdf_params[13] = rotation_inv.xz;
sdf_params[14] = rotation_inv.yx; sdf_params[15] = rotation_inv.yy; sdf_params[16] = rotation_inv.yz;
sdf_params[17] = rotation_inv.zx; sdf_params[18] = rotation_inv.zy; sdf_params[19] = rotation_inv.zz;

// Create kernel and run
Kernel kernel_voxelize_sdf(device, get_N(), "voxelize_sdf", flags, flag, sdf_data, sdf_params);
sdf_data.write_to_device();
sdf_params.write_to_device();
kernel_voxelize_sdf.run();
}
void LBM_Domain::enqueue_unvoxelize_mesh_on_device(const Mesh* mesh, const uchar flag) { // remove voxelized triangle mesh from LBM grid
const float x0=mesh->pmin.x, y0=mesh->pmin.y, z0=mesh->pmin.z, x1=mesh->pmax.x, y1=mesh->pmax.y, z1=mesh->pmax.z; // remove all flags in bounding box of mesh
Kernel kernel_unvoxelize_mesh(device, get_N(), "unvoxelize_mesh", flags, flag, x0, y0, z0, x1, y1, z1);
Expand Down Expand Up @@ -1111,6 +1148,42 @@ void LBM::voxelize_stl(const string& path, const float size, const uchar flag) {
voxelize_stl(path, center(), float3x3(1.0f), size, flag);
}

void LBM::voxelize_sdf(const string& path, const float3& center, const float3x3& rotation, const float scale, const uchar flag) { // voxelize SDF
const SDF* sdf = read_sdf(path);
if(sdf == nullptr) return; // Error already printed by read_sdf

flags.write_to_device();

print_info("Voxelizing SDF on GPU...");
print_info("LBM grid: " + to_string(get_Nx()) + " x " + to_string(get_Ny()) + " x " + to_string(get_Nz()));

// Voxelize on GPU using kernel
for(uint d = 0u; d < get_D(); d++) {
lbm_domain[d]->voxelize_sdf_on_device(sdf, center, rotation, scale, flag);
}
flags.read_from_device();

delete sdf;

// Count total solid voxels
uint total_solid = 0;
for(ulong n = 0; n < get_N(); n++) {
if(flags[n] == flag) total_solid++;
}

flags.write_to_device(); // Push the modified flags back to GPU
print_info("SDF voxelization complete. Total solid voxels: " + to_string(total_solid));
}
void LBM::voxelize_sdf(const string& path, const float3x3& rotation, const float scale, const uchar flag) { // read and voxelize binary SDF file (place in box center)
voxelize_sdf(path, center(), rotation, scale, flag);
}
void LBM::voxelize_sdf(const string& path, const float3& center, const float scale, const uchar flag) { // read and voxelize binary SDF file (no rotation)
voxelize_sdf(path, center, float3x3(1.0f), scale, flag);
}
void LBM::voxelize_sdf(const string& path, const float scale, const uchar flag) { // read and voxelize binary SDF file (place in box center, no rotation)
voxelize_sdf(path, center(), float3x3(1.0f), scale, flag);
}

#ifdef GRAPHICS
int* LBM::Graphics::draw_frame() {
#ifndef UPDATE_FIELDS
Expand Down
5 changes: 5 additions & 0 deletions src/lbm.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -144,6 +144,7 @@ class LBM_Domain {
void set_f(const float fx, const float fy, const float fz) { set_fx(fx); set_fy(fy); set_fz(fz); } // set global froce per volume

void voxelize_mesh_on_device(const Mesh* mesh, const uchar flag=TYPE_S, const float3& rotation_center=float3(0.0f), const float3& linear_velocity=float3(0.0f), const float3& rotational_velocity=float3(0.0f)); // voxelize mesh
void voxelize_sdf_on_device(const SDF* sdf, const float3& center, const float3x3& rotation, const float scale, const uchar flag); // voxelize SDF on GPU
void enqueue_unvoxelize_mesh_on_device(const Mesh* mesh, const uchar flag=TYPE_S); // remove voxelized triangle mesh from LBM grid

#ifdef GRAPHICS
Expand Down Expand Up @@ -542,6 +543,10 @@ class LBM {
void voxelize_stl(const string& path, const float3x3& rotation, const float size=0.0f, const uchar flag=TYPE_S); // read and voxelize binary .stl file (place in box center)
void voxelize_stl(const string& path, const float3& center, const float size=0.0f, const uchar flag=TYPE_S); // read and voxelize binary .stl file (no rotation)
void voxelize_stl(const string& path, const float size=0.0f, const uchar flag=TYPE_S); // read and voxelize binary .stl file (place in box center, no rotation)
void voxelize_sdf(const string& path, const float3& center, const float3x3& rotation, const float scale=1.0f, const uchar flag=TYPE_S); // read and voxelize binary SDF file
void voxelize_sdf(const string& path, const float3x3& rotation, const float scale=1.0f, const uchar flag=TYPE_S); // read and voxelize binary SDF file (place in box center)
void voxelize_sdf(const string& path, const float3& center, const float scale=1.0f, const uchar flag=TYPE_S); // read and voxelize binary SDF file (no rotation)
void voxelize_sdf(const string& path, const float scale=1.0f, const uchar flag=TYPE_S); // read and voxelize binary SDF file (place in box center, no rotation)

#ifdef GRAPHICS
class Graphics {
Expand Down
148 changes: 148 additions & 0 deletions src/utilities.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1211,6 +1211,9 @@ inline float3 operator*(const float3& v, const float3x3& m) { // multiply vector
inline float3 operator*(const float3x3& m, const float3& v) { // multiply matrix with vector
return float3(m.xx*v.x+m.xy*v.y+m.xz*v.z, m.yx*v.x+m.yy*v.y+m.yz*v.z, m.zx*v.x+m.zy*v.y+m.zz*v.z);
}
inline float3x3 transpose(const float3x3& m) { // transpose matrix
return float3x3(m.xx, m.yx, m.zx, m.xy, m.yy, m.zy, m.xz, m.yz, m.zz);
}
inline float3::float3(const float3x3& m) { // extract diagonal of matrix
this->x = m.xx;
this->y = m.yy;
Expand Down Expand Up @@ -4527,6 +4530,151 @@ struct Mesh { // triangle mesh
return fmin(fmin(box_size.x/(pmax.x-pmin.x), box_size.y/(pmax.y-pmin.y)), box_size.z/(pmax.z-pmin.z));
}
};

// ########################### SDF (Signed Distance Field) Support ###########################

struct SDF { // 3D signed distance field
uint Nx = 0u, Ny = 0u, Nz = 0u; // Grid dimensions
float cell_size = 1.0f; // Size of each cell in world units
float3 origin; // World-space origin (lower corner) of the SDF grid
float* data = nullptr; // SDF values: negative=inside, positive=outside, zero=surface

inline SDF(const uint Nx, const uint Ny, const uint Nz, const float cell_size, const float3& origin) {
this->Nx = Nx;
this->Ny = Ny;
this->Nz = Nz;
this->cell_size = cell_size;
this->origin = origin;
this->data = new float[Nx*Ny*Nz];
}

inline ~SDF() {
if(data) delete[] data;
}

inline uint index(const uint x, const uint y, const uint z) const {
return x + Nx * (y + Ny * z);
}

inline float get(const uint x, const uint y, const uint z) const {
return data[index(x, y, z)];
}

inline void set(const uint x, const uint y, const uint z, const float value) {
data[index(x, y, z)] = value;
}

inline float3 get_world_size() const {
return float3((float)Nx * cell_size, (float)Ny * cell_size, (float)Nz * cell_size);
}

inline float3 get_bounding_box_size() const {
return get_world_size();
}

inline float3 get_bounding_box_center() const {
return origin + 0.5f * get_world_size();
}
};

/**
* @brief Read binary SDF file
* @param path Path to binary SDF file
* @return Pointer to SDF structure
*
* File format (binary, little-endian):
* - int32: Nx (grid dimension X)
* - int32: Ny (grid dimension Y)
* - int32: Nz (grid dimension Z)
* - float32: bounds_min_x (bounding box minimum X)
* - float32: bounds_min_y (bounding box minimum Y)
* - float32: bounds_min_z (bounding box minimum Z)
* - float32: bounds_max_x (bounding box maximum X)
* - float32: bounds_max_y (bounding box maximum Y)
* - float32: bounds_max_z (bounding box maximum Z)
* - float32[Nx*Ny*Nz]: SDF values (negative=inside, positive=outside, zero=surface)
* Layout: data[x + Nx*(y + Ny*z)]
*/
inline SDF* read_sdf(const string& path) {
const string filename = create_file_extension(path, ".sdf");
std::ifstream file(filename, std::ios::in | std::ios::binary);

if(!file.is_open()) {
print_error("Error: SDF file \"" + filename + "\" could not be opened!");
return nullptr;
}

// Read header: grid dimensions + bounding box
int Nx, Ny, Nz;
float bounds_min_x, bounds_min_y, bounds_min_z;
float bounds_max_x, bounds_max_y, bounds_max_z;

file.read(reinterpret_cast<char*>(&Nx), sizeof(int));
file.read(reinterpret_cast<char*>(&Ny), sizeof(int));
file.read(reinterpret_cast<char*>(&Nz), sizeof(int));
file.read(reinterpret_cast<char*>(&bounds_min_x), sizeof(float));
file.read(reinterpret_cast<char*>(&bounds_min_y), sizeof(float));
file.read(reinterpret_cast<char*>(&bounds_min_z), sizeof(float));
file.read(reinterpret_cast<char*>(&bounds_max_x), sizeof(float));
file.read(reinterpret_cast<char*>(&bounds_max_y), sizeof(float));
file.read(reinterpret_cast<char*>(&bounds_max_z), sizeof(float));

if(file.fail()) {
print_error("Error: Failed to read SDF header from \"" + filename + "\"!");
file.close();
return nullptr;
}

// Calculate cell size and origin from bounding box
const float3 bounds_min(bounds_min_x, bounds_min_y, bounds_min_z);
const float3 bounds_max(bounds_max_x, bounds_max_y, bounds_max_z);
const float3 bounds_size = bounds_max - bounds_min;

// Cell size = bounding box size / (grid dimensions - 1) for grid corners
// Or = bounding box size / grid dimensions for grid centers
// Using grid centers interpretation:
const float cell_size_x = bounds_size.x / (float)Nx;
const float cell_size_y = bounds_size.y / (float)Ny;
const float cell_size_z = bounds_size.z / (float)Nz;
const float cell_size = (cell_size_x + cell_size_y + cell_size_z) / 3.0f; // Average

const float3 origin = bounds_min;
SDF* sdf = new SDF((uint)Nx, (uint)Ny, (uint)Nz, cell_size, origin);

// Read SDF data
const ulong data_size = (ulong)Nx * (ulong)Ny * (ulong)Nz;
float* temp_data = new float[data_size];
file.read(reinterpret_cast<char*>(temp_data), data_size * sizeof(float));

if(file.fail()) {
print_error("Error: Failed to read SDF data from \"" + filename + "\"!");
delete[] temp_data;
delete sdf;
file.close();
return nullptr;
}

file.close();

// Transpose data from file layout (k varies fastest) to our layout (x varies fastest)
// SDFGen writes: for(i) for(j) for(k) write(value) -> index: i*Ny*Nz + j*Nz + k
// We use: index(x,y,z) = x + Nx*(y + Ny*z) -> x varies fastest
for(uint i = 0; i < (uint)Nx; i++) {
for(uint j = 0; j < (uint)Ny; j++) {
for(uint k = 0; k < (uint)Nz; k++) {
const ulong file_index = (ulong)i * (ulong)Ny * (ulong)Nz + (ulong)j * (ulong)Nz + (ulong)k;
sdf->set(i, j, k, temp_data[file_index]);
}
}
}
delete[] temp_data;

print_info("Loaded SDF: " + to_string(Nx) + "x" + to_string(Ny) + "x" + to_string(Nz) +
" cells, cell_size=" + to_string(cell_size) + "m");

return sdf;
}

inline Mesh* read_stl_raw(const string& path, const bool reposition, const float3& box_size, const float3& center, const float3x3& rotation, const float size) { // read binary .stl file
const string filename = create_file_extension(path, ".stl");
std::ifstream file(filename, std::ios::in|std::ios::binary);
Expand Down