Spaces:
Sleeping
Sleeping
File size: 8,194 Bytes
c83d2d6 | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 | #include <cstdint>
#include <cuda_runtime.h>
struct Point3D {
float x, y, z;
};
struct Triangle {
Point3D v0, v1, v2;
};
struct DistanceSign {
float distance;
int inside;
};
__device__ Point3D cross(const Point3D& v1, const Point3D& v2) {
Point3D result;
result.x = v1.y * v2.z - v1.z * v2.y;
result.y = v1.z * v2.x - v1.x * v2.z;
result.z = v1.x * v2.y - v1.y * v2.x;
return result;
}
__device__ float dot(const Point3D& v1, const Point3D& v2) {
return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}
__device__ Point3D subtract(const Point3D& p1, const Point3D& p2) {
Point3D result = {p1.x - p2.x, p1.y - p2.y, p1.z - p2.z};
return result;
}
__device__ float magnitude(const Point3D &v) {
return sqrtf(v.x * v.x + v.y * v.y + v.z * v.z);
}
__device__ bool is_identical(const Point3D & p1, const Point3D & p2){
Point3D check = subtract(p1, p2);
if(check.x==0 && check.y == 0 && check.z == 0)
return true;
return false;
}
__device__ float squaredDistance(const Point3D& p1, const Point3D& p2) {
return (p1.x - p2.x) * (p1.x - p2.x) +
(p1.y - p2.y) * (p1.y - p2.y) +
(p1.z - p2.z) * (p1.z - p2.z);
}
__device__ Point3D normalize(Point3D v){
float len = sqrtf(dot(v, v));
if (len ==0)
return v;
float scale = 1 / len;
Point3D result = {v.x * scale, v.y * scale, v.z * scale};
return result;
}
__device__ float point_to_line_distance(const Point3D &p, const Point3D &v0, const Point3D &v1) {
// Direction vector of the line
Point3D d = subtract(v1, v0);
// Vector from v0 to point p
Point3D v0_to_p = subtract(p, v0);
// Scalar projection of v0_to_p onto the direction vector d
float t = dot(v0_to_p, d) / dot(d, d);
Point3D closest_point;
// Check where the projection falls
if (t < 0) {
// Projection falls before v0, so the closest point is v0
closest_point = v0;
} else if (t > 1) {
// Projection falls beyond v1, so the closest point is v1
closest_point = v1;
} else {
// Projection falls within the segment, compute the projection point
closest_point.x = v0.x + t * d.x;
closest_point.y = v0.y + t * d.y;
closest_point.z = v0.z + t * d.z;
}
// Calculate the distance between p and the closest point
Point3D closest_to_p = subtract(p, closest_point);
return magnitude(closest_to_p);
}
__device__ DistanceSign pointToTriangleDistanceSign(
const Point3D& queryPoint,
const Point3D& v0,
const Point3D& v1,
const Point3D& v2,
bool inverse=false
) {
Point3D edge0 = subtract(v1, v0);
Point3D edge1 = subtract(v2, v0);
if (is_identical(v0, v1) && is_identical(v0, v2)) {
DistanceSign out = {sqrtf(squaredDistance(queryPoint, v0)), 0};
return out;
}
if (is_identical(v0, v1)) {
DistanceSign out = {point_to_line_distance(queryPoint, v0, v2), 0};
return out;
}
if (is_identical(v0, v2)) {
DistanceSign out = {point_to_line_distance(queryPoint, v0, v1), 0};
return out;
}
if (is_identical(v1, v2)) {
DistanceSign out = {point_to_line_distance(queryPoint, v0, v1), 0};
return out;
}
Point3D normal = cross(edge0, edge1);
if (inverse)
normal = cross(edge1, edge0);
Point3D queryVec = subtract(queryPoint, v0);
if (dot(normal, normal) == 0) {
DistanceSign out = {sqrtf(dot(queryVec, queryVec)), 0};
return out;
}
normal = normalize(normal);
float signedDistanceToPlane = dot(normal, queryVec);
Point3D projectionPoint = {
queryPoint.x - signedDistanceToPlane * normal.x,
queryPoint.y - signedDistanceToPlane * normal.y,
queryPoint.z - signedDistanceToPlane * normal.z
};
edge0 = subtract(v0, v1);
edge1 = subtract(v1, v2);
Point3D edge2 = subtract(v2, v0);
Point3D projVec0 = subtract(v0, projectionPoint);
Point3D projVec1 = subtract(v1, projectionPoint);
Point3D projVec2 = subtract(v2, projectionPoint);
Point3D c0 = cross(edge0, projVec0);
Point3D c1 = cross(edge1, projVec1);
Point3D c2 = cross(edge2, projVec2);
float distance;
if (dot(c0, c1) > 0 && dot(c1, c2) > 0 && dot(c0, c2) > 0) {
distance = fabsf(signedDistanceToPlane);
} else {
distance = 1e6f;
distance = fmin(distance, point_to_line_distance(queryPoint, v0, v1));
distance = fmin(distance, point_to_line_distance(queryPoint, v0, v2));
distance = fmin(distance, point_to_line_distance(queryPoint, v1, v2));
}
DistanceSign out = {distance, signedDistanceToPlane < 0.0f ? 1 : 0};
return out;
}
__device__ void updateSDF(Triangle t, unsigned long long* sdf, const int DIM, const float threshold) {
float minX = fminf(fminf(t.v0.x, t.v1.x), t.v2.x);
float minY = fminf(fminf(t.v0.y, t.v1.y), t.v2.y);
float minZ = fminf(fminf(t.v0.z, t.v1.z), t.v2.z);
float maxX = fmaxf(fmaxf(t.v0.x, t.v1.x), t.v2.x);
float maxY = fmaxf(fmaxf(t.v0.y, t.v1.y), t.v2.y);
float maxZ = fmaxf(fmaxf(t.v0.z, t.v1.z), t.v2.z);
int iMin = max(0, (int)floorf((minX + 0.5) * (DIM-1)));
int jMin = max(0, (int)floorf((minY + 0.5) * (DIM-1)));
int kMin = max(0, (int)floorf((minZ + 0.5) * (DIM-1)));
int iMax = min(DIM - 1, (int)floorf((maxX + 0.5) * (DIM-1)));
int jMax = min(DIM - 1, (int)floorf((maxY + 0.5) * (DIM-1)));
int kMax = min(DIM - 1, (int)floorf((maxZ + 0.5) * (DIM-1)));
int range = (int)(threshold + 1);
iMax = min(DIM - 1, iMax + range);
iMin = max(0, iMin - range);
jMax = min(DIM - 1, jMax + range);
jMin = max(0, jMin - range);
kMax = min(DIM - 1, kMax + range);
kMin = max(0, kMin - range);
float inv_threshold = threshold / DIM;
for (int i = iMin; i <= iMax; ++i) {
for (int j = jMin; j <= jMax; ++j) {
for (int k = kMin; k <= kMax; ++k) {
int idx = i * DIM * DIM + j * DIM + k;
Point3D queryPoint = {
(float)i/(DIM-1) - 0.5f,
(float)j/(DIM-1) - 0.5f,
(float)k/(DIM-1) - 0.5f
};
// Single call: the inverse-normal branch was dead code since
// geometric distance is independent of normal direction.
DistanceSign ds = pointToTriangleDistanceSign(queryPoint, t.v0, t.v1, t.v2);
if (ds.distance < inv_threshold) {
unsigned long long int_dist = (unsigned long long)(ds.distance * 10000000.0f);
// Pack: high bits = distance (drives atomicMin), low bit = inside flag.
// Tie-breaking: inside=0 (outside) wins when distances are equal,
// which is correct for convex edges where both sides are equidistant.
unsigned long long packed = (int_dist << 1) | (unsigned long long)(ds.inside & 1u);
atomicMin(&sdf[idx], packed);
}
}
}
}
}
__global__ void compute_sdf_kernel(float* vertices, int* faces, unsigned long long* sdf, int numTriangles, const int DIM, const float threshold) {
int t = blockIdx.x * blockDim.x + threadIdx.x;
if (t < numTriangles) {
int f0 = faces[t * 3 + 0];
int f1 = faces[t * 3 + 1];
int f2 = faces[t * 3 + 2];
Point3D v0 = {vertices[f0 * 3 + 0], vertices[f0 * 3 + 1], vertices[f0 * 3 + 2]};
Point3D v1 = {vertices[f1 * 3 + 0], vertices[f1 * 3 + 1], vertices[f1 * 3 + 2]};
Point3D v2 = {vertices[f2 * 3 + 0], vertices[f2 * 3 + 1], vertices[f2 * 3 + 2]};
Triangle triangle = {v0, v1, v2};
updateSDF(triangle, sdf, DIM, threshold);
}
}
void compute_valid_sdf_cuda(float* vertices, int* faces, int64_t* sdf, int numTriangles, const int DIM=512, const float threshold=8) {
int blockSize = 256;
int gridSize = (numTriangles + blockSize - 1) / blockSize;
compute_sdf_kernel<<<gridSize, blockSize>>>(
vertices,
faces,
reinterpret_cast<unsigned long long*>(sdf),
numTriangles,
DIM,
threshold
);
}
|