PeriDyno 1.0.0
Loading...
Searching...
No Matches
DistanceField3D.cu
Go to the documentation of this file.
1#include <fstream>
2#include "DistanceField3D.h"
3#include "Vector.h"
4#include "DataTypes.h"
5
6namespace dyno{
7
8 template <typename Coord>
9 __device__ float DistanceToPlane(const Coord &p, const Coord &o, const Coord &n) {
10 return fabs((p - o, n).length());
11 }
12
13 template <typename Coord>
14 __device__ Real DistanceToSegment(Coord& pos, Coord& lo, Coord& hi)
15 {
16 typedef typename Coord::VarType Real;
17 Coord seg = hi - lo;
18 Coord edge1 = pos - lo;
19 Coord edge2 = pos - hi;
20 if (edge1.dot(seg) < 0.0f)
21 {
22 return edge1.norm();
23 }
24 if (edge2.dot(-seg) < 0.0f)
25 {
26 return edge2.norm();
27 }
28 Real length1 = edge1.dot(edge1);
29 seg.normalize();
30 Real length2 = edge1.dot(seg);
31 return std::sqrt(length1 - length2*length2);
32 }
33
34 template <typename Coord>
35 __device__ Real DistanceToSqure(Coord& pos, Coord& lo, Coord& hi, int axis)
36 {
37 typedef typename Coord::VarType Real;
38 Coord n;
39 Coord corner1, corner2, corner3, corner4;
40 Coord loCorner, hiCorner, p;
41 switch (axis)
42 {
43 case 0:
44 corner1 = Coord(lo[0], lo[1], lo[2]);
45 corner2 = Coord(lo[0], hi[1], lo[2]);
46 corner3 = Coord(lo[0], hi[1], hi[2]);
47 corner4 = Coord(lo[0], lo[1], hi[2]);
48 n = Coord(1.0, 0.0, 0.0);
49
50 loCorner = Coord(lo[1], lo[2], 0.0);
51 hiCorner = Coord(hi[1], hi[2], 0.0);
52 p = Coord(pos[1], pos[2], 0.0f);
53 break;
54 case 1:
55 corner1 = Coord(lo[0], lo[1], lo[2]);
56 corner2 = Coord(lo[0], lo[1], hi[2]);
57 corner3 = Coord(hi[0], lo[1], hi[2]);
58 corner4 = Coord(hi[0], lo[1], lo[2]);
59 n = Coord(0.0f, 1.0f, 0.0f);
60
61 loCorner = Coord(lo[0], lo[2], 0.0f);
62 hiCorner = Coord(hi[0], hi[2], 0.0f);
63 p = Coord(pos[0], pos[2], 0.0f);
64 break;
65 case 2:
66 corner1 = Coord(lo[0], lo[1], lo[2]);
67 corner2 = Coord(hi[0], lo[1], lo[2]);
68 corner3 = Coord(hi[0], hi[1], lo[2]);
69 corner4 = Coord(lo[0], hi[1], lo[2]);
70 n = Coord(0.0f, 0.0f, 1.0f);
71
72 loCorner = Coord(lo[0], lo[1], 0.0);
73 hiCorner = Coord(hi[0], hi[1], 0.0);
74 p = Coord(pos[0], pos[1], 0.0f);
75 break;
76 }
77
78 Real dist1 = DistanceToSegment(pos, corner1, corner2);
79 Real dist2 = DistanceToSegment(pos, corner2, corner3);
80 Real dist3 = DistanceToSegment(pos, corner3, corner4);
81 Real dist4 = DistanceToSegment(pos, corner4, corner1);
82 Real dist5 = abs(n.dot(pos - corner1));
83 if (p[0] < hiCorner[0] && p[0] > loCorner[0] && p[1] < hiCorner[1] && p[1] > loCorner[1])
84 return dist5;
85 else
86 return min(min(dist1, dist2), min(dist3, dist4));
87 }
88
89 template <typename Coord>
90 __device__ Real DistanceToBox(Coord& pos, Coord& lo, Coord& hi)
91 {
92 typedef typename Coord::VarType Real;
93 Coord corner0(lo[0], lo[1], lo[2]);
94 Coord corner1(hi[0], lo[1], lo[2]);
95 Coord corner2(hi[0], hi[1], lo[2]);
96 Coord corner3(lo[0], hi[1], lo[2]);
97 Coord corner4(lo[0], lo[1], hi[2]);
98 Coord corner5(hi[0], lo[1], hi[2]);
99 Coord corner6(hi[0], hi[1], hi[2]);
100 Coord corner7(lo[0], hi[1], hi[2]);
101 Real dist0 = (pos - corner0).norm();
102 Real dist1 = (pos - corner1).norm();
103 Real dist2 = (pos - corner2).norm();
104 Real dist3 = (pos - corner3).norm();
105 Real dist4 = (pos - corner4).norm();
106 Real dist5 = (pos - corner5).norm();
107 Real dist6 = (pos - corner6).norm();
108 Real dist7 = (pos - corner7).norm();
109 if (pos[0] < hi[0] && pos[0] > lo[0] && pos[1] < hi[1] && pos[1] > lo[1] && pos[2] < hi[2] && pos[2] > lo[2])
110 {
111 Real distx = min(abs(pos[0] - hi[0]), abs(pos[0] - lo[0]));
112 Real disty = min(abs(pos[1] - hi[1]), abs(pos[1] - lo[1]));
113 Real distz = min(abs(pos[2] - hi[2]), abs(pos[2] - lo[2]));
114 Real mindist = min(distx, disty);
115 mindist = min(mindist, distz);
116 return mindist;
117 }
118 else
119 {
120 Real distx1 = DistanceToSqure(pos, corner0, corner7, 0);
121 Real distx2 = DistanceToSqure(pos, corner1, corner6, 0);
122 Real disty1 = DistanceToSqure(pos, corner0, corner5, 1);
123 Real disty2 = DistanceToSqure(pos, corner3, corner6, 1);
124 Real distz1 = DistanceToSqure(pos, corner0, corner2, 2);
125 Real distz2 = DistanceToSqure(pos, corner4, corner6, 2);
126 return -min(min(min(distx1, distx2), min(disty1, disty2)), min(distz1, distz2));
127 }
128 }
129
130 template <typename Real, typename Coord>
131 __device__ Real DistanceToCylinder(Coord& pos, Coord& center, Real radius, Real height, int axis)
132 {
133 Real distR;
134 Real distH;
135 switch (axis)
136 {
137 case 0:
138 distH = abs(pos[0] - center[0]);
139 distR = Coord(0.0, pos[1] - center[1], pos[2] - center[2]).norm();
140 break;
141 case 1:
142 distH = abs(pos[1] - center[1]);
143 distR = Coord(pos[0] - center[0], 0.0, pos[2] - center[2]).norm();
144 break;
145 case 2:
146 distH = abs(pos[2] - center[2]);
147 distR = Coord(pos[0] - center[0], pos[1] - center[1], 0.0).norm();
148 break;
149 }
150
151 Real halfH = height / 2.0f;
152 if (distH <= halfH && distR <= radius)
153 {
154 return -min(halfH - distH, radius - distR);
155 }
156 else if (distH > halfH && distR <= radius)
157 {
158 return distH - halfH;
159 }
160 else if (distH <= halfH && distR > radius)
161 {
162 return distR - radius;
163 }
164 else
165 {
166// Real l1 = distR - radius;
167// Real l2 = distH - halfH;
168// return sqrt(l1*l1 + l2*l2);
169 return Vector<Real, 2>(distR - radius, distH - halfH).norm();
170 }
171
172
173 }
174
175 template <typename Real, typename Coord>
176 __device__ Real DistanceToSphere(Coord& pos, Coord& center, Real radius)
177 {
178 return (pos - center).length() - radius;
179 }
180
181 template<typename TDataType>
182 DistanceField3D<TDataType>::DistanceField3D()
183 {
184 }
185
186 template<typename TDataType>
187 DistanceField3D<TDataType>::DistanceField3D(std::string filename)
188 {
189 loadSDF(filename);
190 }
191
192 template<typename TDataType>
193 void DistanceField3D<TDataType>::setSpace(const Coord p0, const Coord p1, Real h)
194 {
195 mOrigin = p0;
196 mH = h;
197
198 uint nbx = std::ceil((p1.x - p0.x) / h);
199 uint nby = std::ceil((p1.y - p0.y) / h);
200 uint nbz = std::ceil((p1.z - p0.z) / h);
201
202 mDistances.resize(nbx+1, nby+1, nbz+1);
203 }
204
205 template<typename TDataType>
206 DistanceField3D<TDataType>::~DistanceField3D()
207 {
208 }
209
210 template<typename TDataType>
211 void DistanceField3D<TDataType>::translate(const Coord &t) {
212 mOrigin += t;
213 }
214
215 template <typename Real>
216 __global__ void K_Scale(DArray3D<Real> distance, float s)
217 {
218 int i = threadIdx.x + (blockIdx.x * blockDim.x);
219 int j = threadIdx.y + (blockIdx.y * blockDim.y);
220 int k = threadIdx.z + (blockIdx.z * blockDim.z);
221
222 if (i >= distance.nx()) return;
223 if (j >= distance.ny()) return;
224 if (k >= distance.nz()) return;
225
226 distance(i, j, k) = s*distance(i, j, k);
227 }
228
229 template<typename TDataType>
230 void DistanceField3D<TDataType>::scale(const Real s) {
231 mOrigin[0] *= s;
232 mOrigin[1] *= s;
233 mOrigin[2] *= s;
234 mH *= s;
235
236 dim3 blockSize = make_uint3(8, 8, 8);
237 dim3 gridDims = cudaGridSize3D(make_uint3(mDistances.nx(), mDistances.ny(), mDistances.nz()), blockSize);
238
239 K_Scale << <gridDims, blockSize >> >(mDistances, s);
240 }
241
242 template<typename Real>
243 __global__ void K_Invert(DArray3D<Real> distance)
244 {
245 int i = threadIdx.x + (blockIdx.x * blockDim.x);
246 int j = threadIdx.y + (blockIdx.y * blockDim.y);
247 int k = threadIdx.z + (blockIdx.z * blockDim.z);
248
249 if (i >= distance.nx()) return;
250 if (j >= distance.ny()) return;
251 if (k >= distance.nz()) return;
252
253 distance(i, j, k) = -distance(i, j, k);
254 }
255
256 template<typename TDataType>
257 void DistanceField3D<TDataType>::invertSDF()
258 {
259 dim3 blockSize = make_uint3(8, 8, 8);
260 dim3 gridDims = cudaGridSize3D(make_uint3(mDistances.nx(), mDistances.ny(), mDistances.nz()), blockSize);
261
262 K_Invert << <gridDims, blockSize >> >(mDistances);
263 }
264
265 template <typename Real, typename Coord>
266 __global__ void K_DistanceFieldToBox(DArray3D<Real> distance, Coord start, Real h, Coord lo, Coord hi, bool inverted)
267 {
268 int i = threadIdx.x + (blockIdx.x * blockDim.x);
269 int j = threadIdx.y + (blockIdx.y * blockDim.y);
270 int k = threadIdx.z + (blockIdx.z * blockDim.z);
271
272 if (i >= distance.nx()) return;
273 if (j >= distance.ny()) return;
274 if (k >= distance.nz()) return;
275
276 int sign = inverted ? 1.0f : -1.0f;
277 Coord p = start + Coord(i, j, k) * h;
278
279 distance(i, j, k) = sign*DistanceToBox(p, lo, hi);
280 }
281
282 template<typename TDataType>
283 void DistanceField3D<TDataType>::loadBox(Coord& lo, Coord& hi, bool inverted)
284 {
285 mInverted = inverted;
286
287 dim3 blockSize = make_uint3(4, 4, 4);
288 dim3 gridDims = cudaGridSize3D(make_uint3(mDistances.nx(), mDistances.ny(), mDistances.nz()), blockSize);
289
290 K_DistanceFieldToBox << <gridDims, blockSize >> >(mDistances, mOrigin, mH, lo, hi, inverted);
291 }
292
293 template <typename Real, typename Coord>
294 __global__ void K_DistanceFieldToCylinder(DArray3D<Real> distance, Coord start, Real h, Coord center, Real radius, Real height, int axis, bool inverted)
295 {
296 int i = threadIdx.x + (blockIdx.x * blockDim.x);
297 int j = threadIdx.y + (blockIdx.y * blockDim.y);
298 int k = threadIdx.z + (blockIdx.z * blockDim.z);
299
300 if (i >= distance.nx()) return;
301 if (j >= distance.ny()) return;
302 if (k >= distance.nz()) return;
303
304 int sign = inverted ? -1.0f : 1.0f;
305
306 Coord p = start + Coord(i, j, k)*h;
307
308 distance(i, j, k) = sign*DistanceToCylinder(p, center, radius, height, axis);
309 }
310
311 template<typename TDataType>
312 void DistanceField3D<TDataType>::loadCylinder(Coord& center, Real radius, Real height, int axis, bool inverted)
313 {
314 mInverted = inverted;
315
316 dim3 blockSize = make_uint3(8, 8, 8);
317 dim3 gridDims = cudaGridSize3D(make_uint3(mDistances.nx(), mDistances.ny(), mDistances.nz()), blockSize);
318
319 K_DistanceFieldToCylinder << <gridDims, blockSize >> >(mDistances, mOrigin, mH, center, radius, height, axis, inverted);
320 }
321
322 template <typename Real, typename Coord>
323 __global__ void K_DistanceFieldToSphere(DArray3D<Real> distance, Coord start, Real h, Coord center, Real radius, bool inverted)
324 {
325 int i = threadIdx.x + (blockIdx.x * blockDim.x);
326 int j = threadIdx.y + (blockIdx.y * blockDim.y);
327 int k = threadIdx.z + (blockIdx.z * blockDim.z);
328
329 if (i >= distance.nx()) return;
330 if (j >= distance.ny()) return;
331 if (k >= distance.nz()) return;
332
333 int sign = inverted ? -1.0f : 1.0f;
334
335 Coord p = start + Coord(i, j, k) * h;
336
337 Coord dir = p - center;
338
339 distance(i, j, k) = sign*(dir.norm()-radius);
340 }
341
342 template<typename TDataType>
343 void DistanceField3D<TDataType>::loadSphere(Coord& center, Real radius, bool inverted)
344 {
345 mInverted = inverted;
346
347 dim3 blockSize = make_uint3(8, 8, 8);
348 dim3 gridDims = cudaGridSize3D(make_uint3(mDistances.nx(), mDistances.ny(), mDistances.nz()), blockSize);
349
350 K_DistanceFieldToSphere << <gridDims, blockSize >> >(mDistances, mOrigin, mH, center, radius, inverted);
351 }
352
353 template<typename TDataType>
354 void DistanceField3D<TDataType>::loadSDF(std::string filename, bool inverted)
355 {
356 std::ifstream input(filename.c_str(), std::ios::in);
357 if (!input.is_open())
358 {
359 std::cout << "Reading file " << filename << " error!" << std::endl;
360 return;
361 }
362
363 int nbx, nby, nbz;
364 int xx, yy, zz;
365
366 input >> xx;
367 input >> yy;
368 input >> zz;
369
370 input >> mOrigin[0];
371 input >> mOrigin[1];
372 input >> mOrigin[2];
373
374 Real t_h;
375 input >> t_h;
376
377 std::cout << "SDF: " << xx << ", " << yy << ", " << zz << std::endl;
378 std::cout << "SDF: " << mOrigin[0] << ", " << mOrigin[1] << ", " << mOrigin[2] << std::endl;
379 std::cout << "SDF: " << mOrigin[0] + t_h*xx << ", " << mOrigin[1] + t_h*yy << ", " << mOrigin[2] + t_h*zz << std::endl;
380
381 nbx = xx;
382 nby = yy;
383 nbz = zz;
384 mH = t_h;
385
386 CArray3D<Real> distances(nbx, nby, nbz);
387 for (int k = 0; k < zz; k++) {
388 for (int j = 0; j < yy; j++) {
389 for (int i = 0; i < xx; i++) {
390 float dist;
391 input >> dist;
392 distances(i, j, k) = dist;
393 }
394 }
395 }
396 input.close();
397
398 mDistances.resize(nbx, nby, nbz);
399 mDistances.assign(distances);
400
401 mInverted = inverted;
402 if (inverted)
403 {
404 invertSDF();
405 }
406
407 std::cout << "read data successful" << std::endl;
408 }
409
410 template<typename TDataType>
411 void DistanceField3D<TDataType>::release()
412 {
413 mDistances.clear();
414 }
415
416 template class DistanceField3D<DataType3f>;
417 template class DistanceField3D<DataType3d>;
418}
Generated by doxygen 1.13.2