PeriDyno 1.0.0
Loading...
Searching...
No Matches
PoissonDiskSampling.cu
Go to the documentation of this file.
1#include "GLPointVisualModule.h"
2#include "PoissonDiskSampling.h"
3#include <fstream>
4#include <iostream>
5#include <sstream>
6
7namespace dyno
8{
9 template <typename TDataType>
10 PoissonDiskSampling<TDataType>::PoissonDiskSampling()
11 : SdfSampler<TDataType>()
12 {
13 this->varSpacing()->setRange(0.001, 1.0);
14 area_a = this->varBox_a()->getData() - Coord(this->varSpacing()->getData() * 3);
15 area_b = this->varBox_b()->getData() + Coord(this->varSpacing()->getData() * 3);
16 desired_points = pointNumberRecommend();
17 };
18
19 template <typename TDataType>
20 int PoissonDiskSampling<TDataType>::pointNumberRecommend()
21 {
22 int num = 0;
23 Coord box = area_b - area_a;
24 Real r = this->varSpacing()->getData();
25 return
26 (area_b[0] - area_a[0]) * (area_b[1] - area_a[1]) * (area_b[2] - area_a[2]) / (r * r * r);
27 };
28
29 template <typename TDataType>
30 void PoissonDiskSampling<TDataType>::ConstructGrid()
31 {
32 auto r = this->varSpacing()->getData();
33
34 dx = r / sqrt(3);
35
36
37 if (dx < EPSILON)
38 {
39 std::cout << " The smapling distance in Poisson disk sampling module can not be ZERO!!!! " << std::endl;
40 exit(0);
41 }
42
43 nx = abs(area_b[0] - area_a[0]) / dx;
44 ny = abs(area_b[1] - area_a[1]) / dx;
45 nz = abs(area_b[2] - area_a[2]) / dx;
46
47 gnum = nx * ny * nz;
48
49 m_grid.resize(gnum);
50 for (int i = 0; i < gnum; i++) m_grid[i] = -1;
51 }
52
53
54 template <typename TDataType>
55 GridIndex PoissonDiskSampling<TDataType>::searchGrid(Coord point)
56 {
57 GridIndex index;
58 index.i = (int)((point[0] - area_a[0]) / dx);
59 index.j = (int)((point[1] - area_a[1]) / dx);
60 index.k = (int)((point[2] - area_a[2]) / dx);
61
62 return index;
63 };
64
65 template <typename TDataType>
66 int PoissonDiskSampling<TDataType>::indexTransform(int i, int j, int k)
67 {
68 return i + j * nx + k * nx * ny;
69 };
70
71
72 template <typename TDataType>
73 bool PoissonDiskSampling<TDataType>::collisionJudge2D(Coord point)
74 {
75 bool flag = false;
76 auto r = this->varSpacing()->getData();
77 GridIndex d_index;
78 d_index = searchGrid(point);
79 for (int i = -2; i < 3; i++)
80 for (int j = -2; j < 3; j++)
81 {
82 int mi = d_index.i + i;
83 int mj = d_index.j + j;
84 if ((mi > 0) && (mj > 0) && (mi < nx) && (mj < ny))
85 {
86 int index = indexTransform(mi, mj, d_index.k);
87 if (m_grid[index] != -1)
88 {
89 Coord d = (points[m_grid[index]] - point);
90 if (sqrt(d[0] * d[0] + d[1] * d[1]) - r < EPSILON)
91 {
92 flag = true;
93 }
94 }
95 }
96 }
97 return flag;
98 };
99
100 template <typename TDataType>
101 bool PoissonDiskSampling<TDataType>::collisionJudge(Coord point)
102 {
103
104 bool flag = false;
105 auto r = this->varSpacing()->getData();
106 GridIndex d_index;
107 d_index = searchGrid(point);
108 for (int i = -2; i < 3; i++)
109 for (int j = -2; j < 3; j++)
110 for (int k = -2; k < 3; k++)
111 {
112
113 int mi = d_index.i + i;
114 int mj = d_index.j + j;
115 int mk = d_index.k + k;
116 if ((mi > 0) && (mj > 0) && (mk > 0) && (mi < nx) && (mj < ny) && (mk < nz))
117 {
118 int index = indexTransform(mi, mj, mk);
119 if (m_grid[index] != -1)
120 {
121 Coord d = (points[m_grid[index]] - point);
122 if (sqrt(d[0] * d[0] + d[1] * d[1] + d[2] * d[2]) - r < EPSILON)
123 {
124 flag = true;
125 }
126 }
127 }
128 }
129 return flag;
130 };
131
132
133 template<typename TDataType>
134 std::shared_ptr<DistanceField3D<TDataType>> PoissonDiskSampling<TDataType>::loadSdf()
135 {
136 auto varfilename = this->varSdfFileName()->getData();
137 auto filename = varfilename.string();
138
139 std::shared_ptr<DistanceField3D<TDataType>> sdf = std::make_shared<DistanceField3D<TDataType>>();
140
141 if (filename.size() > 0)
142 {
143
144 if (filename.size() < 5 || filename.substr(filename.size() - 4) != std::string(".sdf")) {
145 std::cerr << "Error: Expected OBJ file with filename of the form <name>.obj.\n";
146 exit(-1);
147 }
148
149 std::ifstream infile(filename);
150 if (!infile) {
151 std::cerr << "Failed to open. Terminating.\n";
152 exit(-1);
153 }
154
155
156
157 sdf->loadSDF(filename, false);
158
159 area_a = sdf->lowerBound();
160 area_b = sdf->upperBound();
161
162 this->varBox_a()->setValue(area_a);
163 this->varBox_b()->setValue(area_b);
164
165 m_h = sdf->getGridSpacing();
166 m_left = sdf->lowerBound();
167
168 std::cout << "SDF is loaded: " << area_a << ", " << area_b << std::endl;
169
170 return sdf;
171 }
172 else
173 {
174 return sdf;
175 }
176 };
177
178 template<typename Real, typename Coord, typename TDataType>
179 __global__ void PDS_PosDetect
180 (
181 DArray<Coord> posArr,
182 DistanceField3D<TDataType> df,
183 DArray<Real> dist
184 )
185 {
186 int pId = threadIdx.x + (blockIdx.x * blockDim.x);
187 if (pId >= posArr.size()) return;
188 Coord pos = posArr[pId];
189 Real temp_dist;
190 Coord normal;
191 df.getDistance(pos, temp_dist, normal);
192 dist[pId] = temp_dist;
193
194 }
195
196
197
198 template <typename TDataType>
199 typename TDataType::Real PoissonDiskSampling<TDataType>::lerp(Real a, Real b, Real alpha)
200 {
201 return (1.0f - alpha) * a + alpha * b;
202 };
203
204
205 template <typename TDataType>
206 typename TDataType::Real PoissonDiskSampling<TDataType>::getDistanceFromSDF(const Coord& p,
207 Coord& normal)
208 {
209 if (!SDF_flag) {
210 return -100;
211 }
212
213 Real d = 0.0f;
214
215 Coord fp = (p - m_left) * Coord(1.0 / m_h, 1.0 / m_h, 1.0 / m_h);
216 const int i = (int)floor(fp[0]);
217 const int j = (int)floor(fp[1]);
218 const int k = (int)floor(fp[2]);
219
220 bool m_bInverted = false;
221
222 if (i < 0 || i >= host_dist.nx() - 1 || j < 0 || j >= host_dist.ny() - 1 || k < 0 || k >= host_dist.nz() - 1) {
223 if (m_bInverted) d = -100000.0f;
224 else d = 100000.0f;
225 normal = Coord(0);
226 return d;
227 }
228
229 Coord ip = Coord(i, j, k);
230
231 Coord alphav = fp - ip;
232 Real alpha = alphav[0];
233 Real beta = alphav[1];
234 Real gamma = alphav[2];
235
236
237 Real d000 = host_dist(i, j, k);
238 Real d100 = host_dist(i + 1, j, k);
239 Real d010 = host_dist(i, j + 1, k);
240 Real d110 = host_dist(i + 1, j + 1, k);
241 Real d001 = host_dist(i, j, k + 1);
242 Real d101 = host_dist(i + 1, j, k + 1);
243 Real d011 = host_dist(i, j + 1, k + 1);
244 Real d111 = host_dist(i + 1, j + 1, k + 1);
245
246 Real dx00 = lerp(d000, d100, alpha);
247 Real dx10 = lerp(d010, d110, alpha);
248 Real dxy0 = lerp(dx00, dx10, beta);
249
250 Real dx01 = lerp(d001, d101, alpha);
251 Real dx11 = lerp(d011, d111, alpha);
252 Real dxy1 = lerp(dx01, dx11, beta);
253
254 Real d0y0 = lerp(d000, d010, beta);
255 Real d0y1 = lerp(d001, d011, beta);
256 Real d0yz = lerp(d0y0, d0y1, gamma);
257
258 Real d1y0 = lerp(d100, d110, beta);
259 Real d1y1 = lerp(d101, d111, beta);
260 Real d1yz = lerp(d1y0, d1y1, gamma);
261
262 Real dx0z = lerp(dx00, dx01, gamma);
263 Real dx1z = lerp(dx10, dx11, gamma);
264
265 normal[0] = d0yz - d1yz;
266 normal[1] = dx0z - dx1z;
267 normal[2] = dxy0 - dxy1;
268
269 Real l = normal.norm();
270 if (l < 0.0001f) normal = Coord(0);
271 else normal = normal.normalize();
272
273 d = (1.0f - gamma) * dxy0 + gamma * dxy1;
274 return d;
275 };
276
277 template<typename TDataType>
278 typename TDataType::Coord PoissonDiskSampling<TDataType>::getOnePointInsideSDF()
279 {
280 Coord normal(0.0f);
281 for (Real ax = area_a[0]; ax < area_b[0]; ax += this->varSpacing()->getData())
282 for (Real ay = area_a[1]; ay < area_b[1]; ay += this->varSpacing()->getData())
283 for (Real az = area_a[2]; az < area_b[2]; az += this->varSpacing()->getData())
284 {
285 Real dr = getDistanceFromSDF(Coord(ax, ay, az), normal);
286 if (dr < 0.0f)
287 {
288 return Coord(ax, ay, az);
289 }
290 }
291 };
292
293
294 template<typename TDataType>
295 void PoissonDiskSampling<TDataType>::resetStates()
296 {
297
298 Real r = this->varSpacing()->getData();
299
300 auto vol = this->getVolume();
301
302 //Adaptive mesh to uniform mesh
303 inputSDF = this->convert2Uniform(vol, r);
304
305 Coord minPoint = inputSDF->lowerBound();
306 Coord maxPoint = inputSDF->upperBound();
307 area_a = minPoint - Coord(r * 3);
308 area_b = maxPoint + Coord(r * 3);
309 this->varBox_a()->setValue(area_a);
310 this->varBox_b()->setValue(area_b);
311 m_h = inputSDF->getGridSpacing();
312 m_left = inputSDF->lowerBound();
313
314 host_dist.resize(inputSDF->nx(), inputSDF->ny(), inputSDF->nz());
315 host_dist.assign(inputSDF->distances());
316
317 SDF_flag = true;
318
319
320 desired_points = pointNumberRecommend();
321
322 std::cout << "Poisson Disk Sampling Start." << std::endl;
323
324 Coord normal(0.0f);
325 this->ConstructGrid();
326 seed_point = (area_a + (area_b - area_a)) / 2;
327
328 /*The seed point must be inside the .SDF boundary */
329 seed_point = getOnePointInsideSDF();
330
331 points.push_back(seed_point);
332 gridIndex = searchGrid(seed_point);
333 int index = indexTransform(gridIndex.i, gridIndex.j, gridIndex.k);
334 m_grid[index] = 0;
335
336 int head = 0;
337 int tail = 1;
338
339 while ((head < desired_points) && (head < tail))
340 {
341 Coord source = points[head];
342 head++;
343 for (int ppp = 0; ppp < attempted_Times; ppp++)
344 {
345 Real theta = (Real)(rand() % 100) / 100.0f;
346 theta = theta * 2 * 3.1415926535;
347
348 Real phi = (Real)(rand() % 100) / 100.0f;
349 phi = phi * 2 * 3.1415926535;
350
351 Real dr = (1 + (Real)(rand() % 100) / 100.0f) * r;
352 Coord offset(0.0f);
353
354 offset = Coord(cos(theta) * sin(phi) * dr, sin(theta) * sin(phi) * dr, cos(phi) * r);
355
356 Coord new_point = source + offset;
357
358 if (((new_point[0] > area_a[0] + r * 3) && (new_point[0] < area_b[0] - r * 3)
359 && (new_point[1] > area_a[1] + r * 3) && (new_point[1] < area_b[1] - r * 3)
360 && (new_point[2] > area_a[2] + r * 3) && (new_point[2] < area_b[2] - r * 3)
361 ))
362 {
363 if (!collisionJudge(new_point) && (tail < desired_points) && (getDistanceFromSDF(new_point, normal) < 0.0f)) {
364 points.push_back(new_point);
365 gridIndex = searchGrid(new_point);
366 m_grid[indexTransform(gridIndex.i, gridIndex.j, gridIndex.k)] = tail;
367 tail++;
368 }
369 }
370 }
371 }
372
373 auto ptSet = this->statePointSet()->getDataPtr();
374 ptSet->setPoints(points);
375 std::cout << "Poisson Disk Sampling is finished." << std::endl;
376
377 points.clear();
378 }
379
380 DEFINE_CLASS(PoissonDiskSampling);
381}
Generated by doxygen 1.13.2