Commit | Line | Data |
---|---|---|
b311480e | 1 | // Created on: 2008-05-30 |
2 | // Created by: Vladislav ROMASHKO | |
3 | // Copyright (c) 2008-2012 OPEN CASCADE SAS | |
4 | // | |
5 | // The content of this file is subject to the Open CASCADE Technology Public | |
6 | // License Version 6.5 (the "License"). You may not use the content of this file | |
7 | // except in compliance with the License. Please obtain a copy of the License | |
8 | // at http://www.opencascade.org and read it completely before using this file. | |
9 | // | |
10 | // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its | |
11 | // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France. | |
12 | // | |
13 | // The Original Code and all software distributed under the License is | |
14 | // distributed on an "AS IS" basis, without warranty of any kind, and the | |
15 | // Initial Developer hereby disclaims all such warranties, including without | |
16 | // limitation, any warranties of merchantability, fitness for a particular | |
17 | // purpose or non-infringement. Please see the License for the specific terms | |
18 | // and conditions governing the rights and limitations under the License. | |
19 | ||
7fd59977 | 20 | |
21 | #include <Voxel_FastConverter.ixx> | |
22 | ||
23 | #include <Bnd_Box.hxx> | |
24 | #include <BRep_Tool.hxx> | |
25 | #include <BRepBndLib.hxx> | |
26 | #include <BRepMesh.hxx> | |
27 | ||
28 | #include <TopoDS.hxx> | |
29 | #include <TopoDS_Face.hxx> | |
30 | #include <TopExp_Explorer.hxx> | |
31 | ||
32 | #include <gp_Lin2d.hxx> | |
33 | #include <gce_MakePln.hxx> | |
34 | ||
35 | #include <ElSLib.hxx> | |
36 | #include <Poly_Triangulation.hxx> | |
37 | #include <IntAna2d_AnaIntersection.hxx> | |
4bee43a9 | 38 | #include <BRepClass3d_SolidClassifier.hxx> |
7fd59977 | 39 | |
40 | // Printing the progress in stdout. | |
41 | //#define CONV_DUMP | |
42 | ||
43 | Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape& shape, | |
44 | Voxel_ROctBoolDS& voxels, | |
45 | const Standard_Real deflection, | |
46 | const Standard_Integer nbx, | |
47 | const Standard_Integer nby, | |
48 | const Standard_Integer nbz, | |
c5e9fb8b P |
49 | const Standard_Integer nbthreads, |
50 | const Standard_Boolean useExistingTriangulation) | |
7fd59977 | 51 | :myShape(shape),myVoxels(&voxels), |
52 | myDeflection(deflection), | |
53 | myNbX(nbx),myNbY(nby),myNbZ(nbz), | |
eafb234b | 54 | myIsBool(2),myNbThreads(nbthreads), |
c5e9fb8b P |
55 | myNbTriangles(0), |
56 | myUseExistingTriangulation(useExistingTriangulation) | |
7fd59977 | 57 | { |
58 | Init(); | |
59 | } | |
60 | ||
61 | Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape& shape, | |
62 | Voxel_BoolDS& voxels, | |
63 | const Standard_Real deflection, | |
64 | const Standard_Integer nbx, | |
65 | const Standard_Integer nby, | |
66 | const Standard_Integer nbz, | |
c5e9fb8b P |
67 | const Standard_Integer nbthreads, |
68 | const Standard_Boolean useExistingTriangulation) | |
7fd59977 | 69 | :myShape(shape),myVoxels(&voxels), |
70 | myDeflection(deflection), | |
71 | myNbX(nbx),myNbY(nby),myNbZ(nbz), | |
eafb234b | 72 | myIsBool(1),myNbThreads(nbthreads), |
c5e9fb8b P |
73 | myNbTriangles(0), |
74 | myUseExistingTriangulation(useExistingTriangulation) | |
7fd59977 | 75 | { |
76 | Init(); | |
77 | } | |
78 | ||
79 | Voxel_FastConverter::Voxel_FastConverter(const TopoDS_Shape& shape, | |
80 | Voxel_ColorDS& voxels, | |
81 | const Standard_Real deflection, | |
82 | const Standard_Integer nbx, | |
83 | const Standard_Integer nby, | |
84 | const Standard_Integer nbz, | |
c5e9fb8b P |
85 | const Standard_Integer nbthreads, |
86 | const Standard_Boolean useExistingTriangulation) | |
7fd59977 | 87 | :myShape(shape),myVoxels(&voxels), |
88 | myDeflection(deflection), | |
89 | myNbX(nbx),myNbY(nby),myNbZ(nbz), | |
eafb234b | 90 | myIsBool(0),myNbThreads(nbthreads), |
c5e9fb8b P |
91 | myNbTriangles(0), |
92 | myUseExistingTriangulation(useExistingTriangulation) | |
7fd59977 | 93 | { |
94 | Init(); | |
95 | } | |
96 | ||
97 | void Voxel_FastConverter::Init() | |
98 | { | |
99 | if (myShape.IsNull()) | |
100 | return; | |
101 | if (myNbThreads < 1) | |
102 | return; | |
103 | ||
104 | // Check number of splits. | |
105 | Voxel_DS* voxels = (Voxel_DS*) myVoxels; | |
106 | if (voxels->GetNbX() != myNbX || voxels->GetNbY() != myNbY || voxels->GetNbZ() != myNbZ) | |
107 | { | |
108 | // Compute boundary box of the shape | |
109 | Bnd_Box box; | |
110 | BRepBndLib::Add(myShape, box); | |
111 | ||
112 | // Define the voxel model by means of the boundary box of shape | |
113 | Standard_Real xmin, ymin, zmin, xmax, ymax, zmax; | |
114 | box.Get(xmin, ymin, zmin, xmax, ymax, zmax); | |
115 | ||
116 | // Initialize the voxels. | |
117 | if (myIsBool == 2) | |
118 | ((Voxel_ROctBoolDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ); | |
119 | else if (myIsBool == 1) | |
120 | ((Voxel_BoolDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ); | |
121 | else if (myIsBool == 0) | |
122 | ((Voxel_ColorDS*) voxels)->Init(xmin, ymin, zmin, xmax - xmin, ymax - ymin, zmax - zmin, myNbX, myNbY, myNbZ); | |
123 | } | |
124 | ||
125 | // Check presence of triangulation. | |
126 | TopLoc_Location L; | |
127 | Standard_Boolean triangulate = Standard_False; | |
128 | TopExp_Explorer expl(myShape, TopAbs_FACE); | |
c5e9fb8b | 129 | if(myUseExistingTriangulation == Standard_False) |
7fd59977 | 130 | { |
c5e9fb8b | 131 | for (; expl.More(); expl.Next()) |
7fd59977 | 132 | { |
c5e9fb8b P |
133 | const TopoDS_Face & F = TopoDS::Face(expl.Current()); |
134 | Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L); | |
135 | if (T.IsNull() || (T->Deflection() > myDeflection)) | |
136 | { | |
137 | triangulate = Standard_True; | |
138 | break; | |
139 | } | |
7fd59977 | 140 | } |
141 | } | |
142 | ||
143 | // Re-create the triangulation. | |
144 | if (triangulate) | |
145 | { | |
146 | BRepMesh::Mesh(myShape, myDeflection); | |
147 | } | |
148 | ||
149 | // Compute the number of triangles. | |
150 | myNbTriangles = 0; | |
151 | expl.Init(myShape, TopAbs_FACE); | |
152 | for (; expl.More(); expl.Next()) | |
153 | { | |
034b4775 | 154 | const TopoDS_Face & F = TopoDS::Face(expl.Current()); |
7fd59977 | 155 | Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L); |
034b4775 P |
156 | if (T.IsNull() == Standard_False) |
157 | myNbTriangles += T->NbTriangles(); | |
7fd59977 | 158 | } |
159 | } | |
160 | ||
161 | // Destructor | |
162 | void Voxel_FastConverter::Destroy() | |
163 | { | |
164 | ||
165 | } | |
166 | ||
167 | Standard_Boolean Voxel_FastConverter::Convert(Standard_Integer& progress, | |
168 | const Standard_Integer ithread) | |
169 | { | |
170 | if (ithread == 1) | |
171 | progress = 0; | |
172 | #ifdef CONV_DUMP | |
173 | if (ithread == 1) | |
174 | printf("Progress = %d \r", progress); | |
175 | #endif | |
176 | ||
177 | if (myNbX <= 0 || myNbY <= 0 || myNbZ <= 0) | |
178 | return Standard_False; | |
179 | ||
c5e9fb8b P |
180 | if(myNbTriangles == 0) |
181 | return Standard_False; | |
182 | ||
7fd59977 | 183 | // Half of diagonal of a voxel |
184 | Voxel_DS* ds = (Voxel_DS*) myVoxels; | |
185 | Standard_Real dx = ds->GetXLen() / (Standard_Real) ds->GetNbX(), | |
186 | dy = ds->GetYLen() / (Standard_Real) ds->GetNbY(), | |
187 | dz = ds->GetZLen() / (Standard_Real) ds->GetNbZ(); | |
188 | Standard_Real hdiagonal = sqrt(dx * dx + dy * dy + dz * dz); | |
189 | hdiagonal /= 2.0; | |
190 | ||
191 | // Compute the scope of triangles for current thread | |
192 | Standard_Integer start_thread_triangle = 1, end_thread_triangle = myNbTriangles, ithread_triangle = 0; | |
03679c48 P |
193 | if(myNbTriangles < myNbThreads) |
194 | { | |
195 | if(ithread != 1) | |
196 | return Standard_False; | |
197 | //in case we're in thread one process all triangles | |
198 | } | |
199 | else | |
200 | { | |
201 | div_t division = div(myNbTriangles, myNbThreads); | |
202 | start_thread_triangle = (ithread - 1) * division.quot + 1; | |
203 | end_thread_triangle = (ithread - 0) * division.quot; | |
204 | ||
205 | if(ithread == myNbThreads) | |
206 | end_thread_triangle += division.rem; | |
207 | } | |
7fd59977 | 208 | |
209 | // Convert | |
210 | TopLoc_Location L; | |
302f96fb | 211 | Standard_Integer iprogress = 0; |
7fd59977 | 212 | Standard_Integer n1, n2, n3; |
213 | Standard_Integer ixmin, iymin, izmin, ixmax, iymax, izmax; | |
214 | Standard_Real xmin, ymin, zmin, xmax, ymax, zmax; | |
215 | TopExp_Explorer expl(myShape, TopAbs_FACE); | |
216 | for (; expl.More(); expl.Next()) | |
217 | { | |
034b4775 | 218 | const TopoDS_Face & F = TopoDS::Face(expl.Current()); |
7fd59977 | 219 | Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L); |
034b4775 P |
220 | if (T.IsNull()) |
221 | continue; | |
7fd59977 | 222 | |
223 | gp_Trsf trsf; | |
224 | Standard_Boolean transform = !L.IsIdentity(); | |
225 | if (transform) | |
226 | trsf = L.Transformation(); | |
227 | ||
228 | const TColgp_Array1OfPnt& nodes = T->Nodes(); | |
229 | const Poly_Array1OfTriangle& triangles = T->Triangles(); | |
230 | Standard_Integer itriangle = triangles.Lower(), nb_triangles = triangles.Upper(); | |
231 | for (; itriangle <= nb_triangles; itriangle++) | |
232 | { | |
233 | ithread_triangle++; | |
da8536ad P |
234 | if (ithread_triangle < start_thread_triangle ) |
235 | continue; | |
236 | if (ithread_triangle > end_thread_triangle) | |
237 | { | |
238 | if (ithread == 1) | |
239 | progress = 100; | |
240 | #ifdef CONV_DUMP | |
241 | if (ithread == 1) | |
242 | printf("Progress = %d \r", progress); | |
243 | #endif | |
244 | return Standard_True; | |
245 | } | |
7fd59977 | 246 | |
247 | const Poly_Triangle& t = triangles.Value(itriangle); | |
248 | t.Get(n1, n2, n3); | |
249 | gp_Pnt p1 = nodes.Value(n1); | |
250 | gp_Pnt p2 = nodes.Value(n2); | |
251 | gp_Pnt p3 = nodes.Value(n3); | |
252 | if (transform) | |
253 | { | |
da8536ad P |
254 | p1.Transform(trsf); |
255 | p2.Transform(trsf); | |
256 | p3.Transform(trsf); | |
7fd59977 | 257 | } |
258 | ||
259 | // Get boundary box of the triangle | |
260 | GetBndBox(p1, p2, p3, xmin, ymin, zmin, xmax, ymax, zmax); | |
261 | ||
262 | // Find the range of voxels inside the boudary box of the triangle. | |
263 | if (!ds->GetVoxel(xmin, ymin, zmin, ixmin, iymin, izmin)) | |
264 | continue; | |
265 | if (!ds->GetVoxel(xmax, ymax, zmax, ixmax, iymax, izmax)) | |
266 | continue; | |
267 | ||
268 | // Refuse voxels for whom distance from their center to plane of triangle is greater than half of diagonal. | |
269 | // Make a line from center of each voxel to the center of triangle and | |
270 | // compute intersection of the line with sides of triangle. | |
271 | // Refuse the voxel in case of intersection. | |
272 | gce_MakePln mkPlane(p1, p2, p3); | |
273 | if (!mkPlane.IsDone()) | |
274 | continue; | |
275 | gp_Pln plane = mkPlane.Value(); | |
276 | ComputeVoxelsNearTriangle(plane, p1, p2, p3, hdiagonal, ixmin, iymin, izmin, ixmax, iymax, izmax); | |
277 | ||
278 | // Progress | |
279 | if (ithread == 1) | |
280 | { | |
281 | iprogress++; | |
282 | progress = (Standard_Integer) ( (Standard_Real) iprogress / (Standard_Real) myNbTriangles * 100.0 ); | |
283 | } | |
284 | #ifdef CONV_DUMP | |
285 | if (ithread == 1 && prev_progress != progress) | |
286 | { | |
287 | printf("Progress = %d \r", progress); | |
288 | prev_progress = progress; | |
289 | } | |
290 | #endif | |
291 | ||
292 | } // iteration of triangles | |
293 | } // iteration of faces | |
294 | ||
295 | if (ithread == 1) | |
296 | progress = 100; | |
297 | #ifdef CONV_DUMP | |
298 | if (ithread == 1) | |
299 | printf("Progress = %d \r", progress); | |
300 | #endif | |
301 | return Standard_True; | |
302 | } | |
303 | ||
f03671a0 P |
304 | Standard_Boolean Voxel_FastConverter::ConvertUsingSAT(Standard_Integer& progress, |
305 | const Standard_Integer ithread) | |
306 | { | |
307 | if (ithread == 1) | |
308 | progress = 0; | |
309 | #ifdef CONV_DUMP | |
310 | if (ithread == 1) | |
311 | printf("Progress = %d \r", progress); | |
312 | #endif | |
313 | ||
314 | if (myNbX <= 0 || myNbY <= 0 || myNbZ <= 0) | |
315 | return Standard_False; | |
316 | ||
317 | if(myNbTriangles == 0) | |
318 | return Standard_False; | |
319 | ||
320 | // Half of size of a voxel (also for Voxel_ROctBoolDS) | |
321 | Voxel_DS* ds = (Voxel_DS*) myVoxels; | |
322 | Standard_Real dx = ds->GetXLen() / (Standard_Real) ds->GetNbX(), | |
323 | dy = ds->GetYLen() / (Standard_Real) ds->GetNbY(), | |
324 | dz = ds->GetZLen() / (Standard_Real) ds->GetNbZ(); | |
325 | gp_Pnt extents(dx/2.0, dy/2.0, dz/2.0); | |
326 | gp_Pnt extents2(dx/4.0, dy/4.0, dz/4.0); | |
327 | gp_Pnt extents4(dx/8.0, dy/8.0, dz/8.0); | |
328 | ||
329 | // Compute the scope of triangles for current thread | |
330 | Standard_Integer start_thread_triangle = 1, end_thread_triangle = myNbTriangles, ithread_triangle = 0; | |
03679c48 P |
331 | if(myNbTriangles < myNbThreads) |
332 | { | |
333 | if(ithread != 1) | |
334 | return Standard_False; | |
335 | //in case we're in thread one process all triangles | |
336 | } | |
337 | else | |
338 | { | |
339 | div_t division = div(myNbTriangles, myNbThreads); | |
340 | start_thread_triangle = (ithread - 1) * division.quot + 1; | |
341 | end_thread_triangle = (ithread - 0) * division.quot; | |
342 | ||
343 | if(ithread == myNbThreads) | |
344 | end_thread_triangle += division.rem; | |
345 | } | |
f03671a0 P |
346 | |
347 | // Convert | |
348 | TopLoc_Location L; | |
302f96fb | 349 | Standard_Integer iprogress = 0; |
f03671a0 P |
350 | Standard_Integer n1, n2, n3; |
351 | Standard_Integer ixmin, iymin, izmin, ixmax, iymax, izmax; | |
352 | Standard_Real xmin, ymin, zmin, xmax, ymax, zmax; | |
353 | TopExp_Explorer expl(myShape, TopAbs_FACE); | |
354 | for (; expl.More(); expl.Next()) | |
355 | { | |
356 | const TopoDS_Face & F = TopoDS::Face(expl.Current()); | |
357 | Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L); | |
358 | if (T.IsNull()) | |
359 | continue; | |
360 | ||
361 | gp_Trsf trsf; | |
362 | Standard_Boolean transform = !L.IsIdentity(); | |
363 | if (transform) | |
364 | trsf = L.Transformation(); | |
365 | ||
366 | const TColgp_Array1OfPnt& nodes = T->Nodes(); | |
367 | const Poly_Array1OfTriangle& triangles = T->Triangles(); | |
368 | Standard_Integer itriangle = triangles.Lower(), nb_triangles = triangles.Upper(); | |
369 | for (; itriangle <= nb_triangles; itriangle++) | |
370 | { | |
371 | ithread_triangle++; | |
372 | if (ithread_triangle < start_thread_triangle ) | |
373 | continue; | |
374 | if (ithread_triangle > end_thread_triangle) | |
375 | { | |
376 | if (ithread == 1) | |
377 | progress = 100; | |
378 | #ifdef CONV_DUMP | |
379 | if (ithread == 1) | |
380 | printf("Progress = %d \r", progress); | |
381 | #endif | |
382 | return Standard_True; | |
383 | } | |
384 | ||
385 | const Poly_Triangle& t = triangles.Value(itriangle); | |
386 | ||
387 | t.Get(n1, n2, n3); | |
388 | gp_Pnt p1 = nodes.Value(n1); | |
389 | gp_Pnt p2 = nodes.Value(n2); | |
390 | gp_Pnt p3 = nodes.Value(n3); | |
391 | if (transform) | |
392 | { | |
393 | p1.Transform(trsf); | |
394 | p2.Transform(trsf); | |
395 | p3.Transform(trsf); | |
396 | } | |
397 | ||
398 | // Get boundary box of the triangle | |
399 | GetBndBox(p1, p2, p3, xmin, ymin, zmin, xmax, ymax, zmax); | |
400 | ||
401 | // Find the range of voxels inside the boudary box of the triangle. | |
402 | if (!ds->GetVoxel(xmin, ymin, zmin, ixmin, iymin, izmin)) | |
403 | continue; | |
404 | if (!ds->GetVoxel(xmax, ymax, zmax, ixmax, iymax, izmax)) | |
405 | continue; | |
406 | ||
407 | // Perform triangle-box intersection to find the voxels resulting from the processed triangle.; | |
408 | // Using SAT theorem to quickly find the intersection. | |
409 | ComputeVoxelsNearTriangle(p1, p2, p3, | |
410 | extents, extents2, extents4, | |
411 | ixmin, iymin, izmin, ixmax, iymax, izmax); | |
412 | ||
413 | // Progress | |
414 | if (ithread == 1) | |
415 | { | |
416 | iprogress++; | |
417 | progress = (Standard_Integer) ( (Standard_Real) iprogress / (Standard_Real) myNbTriangles * 100.0 ); | |
418 | } | |
419 | #ifdef CONV_DUMP | |
420 | if (ithread == 1 && prev_progress != progress) | |
421 | { | |
422 | printf("Progress = %d \r", progress); | |
423 | prev_progress = progress; | |
424 | } | |
425 | #endif | |
426 | ||
427 | } // iteration of triangles | |
428 | } // iteration of faces | |
429 | ||
430 | if (ithread == 1) | |
431 | progress = 100; | |
432 | #ifdef CONV_DUMP | |
433 | if (ithread == 1) | |
434 | printf("Progress = %d \r", progress); | |
435 | #endif | |
436 | return Standard_True; | |
437 | } | |
438 | ||
7fd59977 | 439 | Standard_Boolean Voxel_FastConverter::FillInVolume(const Standard_Byte inner, |
35e08fe8 | 440 | const Standard_Integer /*ithread*/) |
7fd59977 | 441 | { |
442 | Voxel_DS* ds = (Voxel_DS*) myVoxels; | |
443 | Standard_Integer ix, iy, iz, nbx = ds->GetNbX(), nby = ds->GetNbY(), nbz = ds->GetNbZ(); | |
444 | Standard_Boolean prev_surface, surface, volume; | |
445 | ||
446 | if (inner) | |
447 | { | |
448 | // Fill-in internal voxels by the value "inner" | |
449 | for (ix = 0; ix < nbx; ix++) | |
450 | { | |
451 | for (iy = 0; iy < nby; iy++) | |
452 | { | |
453 | // Check existence of volume. | |
454 | volume = Standard_False; | |
455 | surface = Standard_False; | |
456 | prev_surface = Standard_False; | |
457 | for (iz = 0; iz < nbz; iz++) | |
458 | { | |
459 | surface = (myIsBool == 1) ? | |
460 | ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True : | |
461 | ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0; | |
462 | if (prev_surface && !surface) | |
463 | { | |
464 | volume = !volume; | |
465 | } | |
466 | prev_surface = surface; | |
467 | } | |
468 | if (volume) | |
469 | continue; | |
470 | ||
471 | // Fill-in the volume. | |
472 | volume = Standard_False; | |
473 | surface = Standard_False; | |
474 | prev_surface = Standard_False; | |
475 | for (iz = 0; iz < nbz; iz++) | |
476 | { | |
477 | surface = (myIsBool == 1) ? | |
478 | ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True : | |
479 | ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0; | |
480 | if (prev_surface && !surface) | |
481 | { | |
482 | volume = !volume; | |
483 | } | |
484 | if (volume && !surface) | |
485 | { | |
486 | (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) : | |
487 | ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner); | |
488 | } | |
489 | prev_surface = surface; | |
490 | } | |
491 | } | |
492 | } | |
493 | } | |
494 | else | |
495 | { | |
496 | // Set value of interbal voxels to 0 ("inner" = 0) | |
497 | Standard_Boolean next_surface; | |
498 | for (ix = 0; ix < nbx; ix++) | |
499 | { | |
500 | for (iy = 0; iy < nby; iy++) | |
501 | { | |
502 | volume = Standard_False; | |
503 | surface = Standard_False; | |
504 | prev_surface = Standard_False; | |
505 | next_surface = Standard_False; | |
506 | for (iz = 0; iz < nbz; iz++) | |
507 | { | |
508 | surface = (myIsBool == 1) ? | |
509 | ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True : | |
510 | ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0; | |
511 | if (prev_surface != surface) | |
512 | { | |
513 | volume = !volume; | |
514 | } | |
515 | if (volume && iz + 1 < nbz) | |
516 | { | |
517 | next_surface = (myIsBool == 1) ? | |
518 | ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz + 1) == Standard_True : | |
519 | ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz + 1) > 0; | |
520 | } | |
521 | if (volume && prev_surface == surface && next_surface) | |
522 | { | |
523 | (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) : | |
524 | ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner); | |
525 | } | |
526 | prev_surface = surface; | |
527 | } | |
528 | } | |
529 | } | |
530 | } | |
531 | ||
532 | return Standard_True; | |
533 | } | |
534 | ||
35e08fe8 | 535 | Standard_Boolean Voxel_FastConverter::FillInVolume(const Standard_Byte inner, const TopoDS_Shape & shape, const Standard_Integer /*ithread*/) |
4bee43a9 P |
536 | { |
537 | Voxel_DS* ds = (Voxel_DS*) myVoxels; | |
538 | Standard_Integer ix, iy, iz, nbx = ds->GetNbX(), nby = ds->GetNbY(), nbz = ds->GetNbZ(); | |
539 | Standard_Boolean prev_surface, surface, volume, isOnVerticalSurface; | |
540 | ||
541 | BRepClass3d_SolidClassifier solidClassifier(shape); | |
542 | Standard_Real xc, yc, zc; | |
543 | ||
544 | if (inner) | |
545 | { | |
546 | // Fill-in internal voxels by the value "inner" | |
547 | for (ix = 0; ix < nbx; ix++) | |
548 | { | |
549 | for (iy = 0; iy < nby; iy++) | |
550 | { | |
551 | // Check existence of volume. | |
552 | volume = Standard_False; | |
553 | surface = Standard_False; | |
554 | prev_surface = Standard_False; | |
555 | isOnVerticalSurface = Standard_False; | |
556 | for (iz = 0; iz < nbz; iz++) | |
557 | { | |
558 | surface = (myIsBool == 1) ? | |
559 | ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True : | |
560 | ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0; | |
561 | if (prev_surface && !surface) | |
562 | { | |
563 | if(isOnVerticalSurface) | |
564 | { | |
565 | isOnVerticalSurface = Standard_False; | |
566 | ((Voxel_BoolDS*)myVoxels)->GetCenter(ix, iy, iz, xc, yc, zc); | |
567 | gp_Pnt P(xc, yc, zc); | |
568 | solidClassifier.Perform(P, Precision::Confusion()); | |
569 | ||
570 | if(solidClassifier.State() == TopAbs_IN) | |
571 | volume = Standard_True; | |
572 | else | |
573 | volume = Standard_False; | |
574 | } | |
575 | else | |
576 | volume = !volume; | |
577 | } | |
578 | if(prev_surface && surface) | |
579 | isOnVerticalSurface = Standard_True; | |
580 | else | |
581 | isOnVerticalSurface = Standard_False; | |
582 | prev_surface = surface; | |
583 | } | |
584 | if (volume) | |
585 | continue; | |
586 | ||
587 | // Fill-in the volume. | |
588 | volume = Standard_False; | |
589 | surface = Standard_False; | |
590 | prev_surface = Standard_False; | |
591 | isOnVerticalSurface = Standard_False; | |
592 | for (iz = 0; iz < nbz; iz++) | |
593 | { | |
594 | surface = (myIsBool == 1) ? | |
595 | ((Voxel_BoolDS*)myVoxels)->Get(ix, iy, iz) == Standard_True : | |
596 | ((Voxel_ColorDS*)myVoxels)->Get(ix, iy, iz) > 0; | |
597 | if (prev_surface && !surface) | |
598 | { | |
599 | if(isOnVerticalSurface) | |
600 | { | |
601 | isOnVerticalSurface = Standard_False; | |
602 | ((Voxel_BoolDS*)myVoxels)->GetCenter(ix, iy, iz, xc, yc, zc); | |
603 | gp_Pnt P(xc, yc, zc); | |
604 | solidClassifier.Perform(P, Precision::Confusion()); | |
605 | ||
606 | if(solidClassifier.State() == TopAbs_IN) | |
607 | volume = Standard_True; | |
608 | else | |
609 | volume = Standard_False; | |
610 | } | |
611 | else | |
612 | volume = !volume; | |
613 | } | |
614 | if (volume && !surface) | |
615 | { | |
616 | (myIsBool == 1) ? ((Voxel_BoolDS*)myVoxels)->Set(ix, iy, iz, inner) : | |
617 | ((Voxel_ColorDS*)myVoxels)->Set(ix, iy, iz, inner); | |
618 | } | |
619 | if(prev_surface && surface) | |
620 | isOnVerticalSurface = Standard_True; | |
621 | else | |
622 | isOnVerticalSurface = Standard_False; | |
623 | prev_surface = surface; | |
624 | } | |
625 | } | |
626 | } | |
627 | } | |
628 | ||
629 | return Standard_True; | |
630 | } | |
631 | ||
7fd59977 | 632 | void Voxel_FastConverter::GetBndBox(const gp_Pnt& p1, |
633 | const gp_Pnt& p2, | |
634 | const gp_Pnt& p3, | |
635 | Standard_Real& xmin, | |
636 | Standard_Real& ymin, | |
637 | Standard_Real& zmin, | |
638 | Standard_Real& xmax, | |
639 | Standard_Real& ymax, | |
640 | Standard_Real& zmax) const | |
641 | { | |
642 | // P1: | |
643 | xmin = p1.X(); | |
644 | ymin = p1.Y(); | |
645 | zmin = p1.Z(); | |
646 | xmax = p1.X(); | |
647 | ymax = p1.Y(); | |
648 | zmax = p1.Z(); | |
649 | // P2: | |
650 | if (xmin > p2.X()) | |
651 | xmin = p2.X(); | |
652 | if (ymin > p2.Y()) | |
653 | ymin = p2.Y(); | |
654 | if (zmin > p2.Z()) | |
655 | zmin = p2.Z(); | |
656 | if (xmax < p2.X()) | |
657 | xmax = p2.X(); | |
658 | if (ymax < p2.Y()) | |
659 | ymax = p2.Y(); | |
660 | if (zmax < p2.Z()) | |
661 | zmax = p2.Z(); | |
662 | // P3: | |
663 | if (xmin > p3.X()) | |
664 | xmin = p3.X(); | |
665 | if (ymin > p3.Y()) | |
666 | ymin = p3.Y(); | |
667 | if (zmin > p3.Z()) | |
668 | zmin = p3.Z(); | |
669 | if (xmax < p3.X()) | |
670 | xmax = p3.X(); | |
671 | if (ymax < p3.Y()) | |
672 | ymax = p3.Y(); | |
673 | if (zmax < p3.Z()) | |
674 | zmax = p3.Z(); | |
675 | } | |
676 | ||
677 | // This method is copied from Voxel_ShapeIntersector.cxx | |
678 | static Standard_Boolean mayIntersect(const gp_Pnt2d& p11, const gp_Pnt2d& p12, | |
679 | const gp_Pnt2d& p21, const gp_Pnt2d& p22) | |
680 | { | |
681 | if (p11.X() > p21.X() && p11.X() > p22.X() && p12.X() > p21.X() && p12.X() > p22.X()) | |
682 | return Standard_False; | |
683 | if (p11.X() < p21.X() && p11.X() < p22.X() && p12.X() < p21.X() && p12.X() < p22.X()) | |
684 | return Standard_False; | |
685 | if (p11.Y() > p21.Y() && p11.Y() > p22.Y() && p12.Y() > p21.Y() && p12.Y() > p22.Y()) | |
686 | return Standard_False; | |
687 | if (p11.Y() < p21.Y() && p11.Y() < p22.Y() && p12.Y() < p21.Y() && p12.Y() < p22.Y()) | |
688 | return Standard_False; | |
689 | return Standard_True; | |
690 | } | |
691 | ||
692 | void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pln& plane, | |
693 | const gp_Pnt& p1, | |
694 | const gp_Pnt& p2, | |
695 | const gp_Pnt& p3, | |
696 | const Standard_Real hdiagonal, | |
697 | const Standard_Integer ixmin, | |
698 | const Standard_Integer iymin, | |
699 | const Standard_Integer izmin, | |
700 | const Standard_Integer ixmax, | |
701 | const Standard_Integer iymax, | |
702 | const Standard_Integer izmax) const | |
703 | { | |
704 | gp_Pnt pc; | |
705 | Standard_Real xc, yc, zc, uc, vc, u1, v1, u2, v2, u3, v3; | |
706 | Standard_Integer ix, iy, iz; | |
707 | IntAna2d_AnaIntersection intersector2d; | |
708 | ||
709 | // Project points of triangle onto the plane | |
710 | ElSLib::Parameters(plane, p1, u1, v1); | |
711 | ElSLib::Parameters(plane, p2, u2, v2); | |
712 | ElSLib::Parameters(plane, p3, u3, v3); | |
713 | ||
714 | // Make lines of triangle | |
715 | gp_Pnt2d p2d1(u1, v1), p2d2(u2, v2), p2d3(u3, v3), p2dt((u1+u2+u3)/3.0,(v1+v2+v3)/3.0), p2dc; | |
716 | gp_Vec2d v2d12(p2d1, p2d2), v2d23(p2d2, p2d3), v2d31(p2d3, p2d1); | |
717 | gp_Lin2d L1(p2d1, v2d12), L2(p2d2, v2d23), L3(p2d3, v2d31), Lv; | |
718 | Standard_Real d1 = p2d1.Distance(p2d2) - Precision::Confusion(), | |
719 | d2 = p2d2.Distance(p2d3) - Precision::Confusion(), | |
720 | d3 = p2d3.Distance(p2d1) - Precision::Confusion(), dv; | |
721 | ||
722 | Voxel_DS* ds = (Voxel_DS*) myVoxels; | |
723 | for (ix = ixmin; ix <= ixmax; ix++) | |
724 | { | |
725 | for (iy = iymin; iy <= iymax; iy++) | |
726 | { | |
727 | for (iz = izmin; iz <= izmax; iz++) | |
728 | { | |
729 | ds->GetCenter(ix, iy, iz, xc, yc, zc); | |
730 | pc.SetCoord(xc, yc, zc); | |
731 | if (plane.Distance(pc) < hdiagonal) | |
732 | { | |
733 | ElSLib::Parameters(plane, pc, uc, vc); | |
734 | p2dc.SetCoord(uc, vc); | |
735 | ||
736 | gp_Vec2d v2dct(p2dc, p2dt); | |
737 | dv = v2dct.Magnitude() - Precision::Confusion(); | |
738 | Lv.SetLocation(p2dc); | |
739 | Lv.SetDirection(v2dct); | |
740 | ||
741 | // Side 1: | |
742 | if (mayIntersect(p2d1, p2d2, p2dc, p2dt)) | |
743 | { | |
744 | intersector2d.Perform(Lv, L1); | |
745 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
746 | { | |
747 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
748 | Standard_Real param1 = i2d.ParamOnFirst(); | |
749 | Standard_Real param2 = i2d.ParamOnSecond(); | |
750 | if (param1 > Precision::Confusion() && param1 < dv && | |
751 | param2 > Precision::Confusion() && param2 < d1) | |
752 | { | |
753 | continue; | |
754 | } | |
755 | } | |
756 | } | |
757 | ||
758 | // Side 2: | |
759 | if (mayIntersect(p2d2, p2d3, p2dc, p2dt)) | |
760 | { | |
761 | intersector2d.Perform(Lv, L2); | |
762 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
763 | { | |
764 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
765 | Standard_Real param1 = i2d.ParamOnFirst(); | |
766 | Standard_Real param2 = i2d.ParamOnSecond(); | |
767 | if (param1 > Precision::Confusion() && param1 < dv && | |
768 | param2 > Precision::Confusion() && param2 < d2) | |
769 | { | |
770 | continue; | |
771 | } | |
772 | } | |
773 | } | |
774 | ||
775 | // Side 3: | |
776 | if (mayIntersect(p2d3, p2d1, p2dc, p2dt)) | |
777 | { | |
778 | intersector2d.Perform(Lv, L3); | |
779 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
780 | { | |
781 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
782 | Standard_Real param1 = i2d.ParamOnFirst(); | |
783 | Standard_Real param2 = i2d.ParamOnSecond(); | |
784 | if (param1 > Precision::Confusion() && param1 < dv && | |
785 | param2 > Precision::Confusion() && param2 < d3) | |
786 | { | |
787 | continue; | |
788 | } | |
789 | } | |
790 | } | |
791 | ||
792 | // Set positive value to this voxel: | |
793 | switch (myIsBool) | |
794 | { | |
795 | case 0: | |
796 | ((Voxel_ColorDS*) myVoxels)->Set(ix, iy, iz, 15); | |
797 | break; | |
798 | case 1: | |
799 | ((Voxel_BoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True); | |
800 | break; | |
801 | case 2: | |
802 | { | |
803 | //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True); | |
804 | ||
805 | // Check intersection between the triangle & sub-voxels of the voxel. | |
806 | Standard_Real hdiagonal2 = hdiagonal / 2.0, hdiagonal4 = hdiagonal / 4.0; | |
807 | for (Standard_Integer i = 0; i < 8; i++) | |
808 | { | |
809 | ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, xc, yc, zc); | |
810 | pc.SetCoord(xc, yc, zc); | |
811 | if (plane.Distance(pc) < hdiagonal2) | |
812 | { | |
813 | ElSLib::Parameters(plane, pc, uc, vc); | |
814 | p2dc.SetCoord(uc, vc); | |
815 | ||
816 | gp_Vec2d v2dct(p2dc, p2dt); | |
817 | dv = v2dct.Magnitude() - Precision::Confusion(); | |
818 | Lv.SetLocation(p2dc); | |
819 | Lv.SetDirection(v2dct); | |
820 | ||
821 | // Side 1: | |
822 | if (mayIntersect(p2d1, p2d2, p2dc, p2dt)) | |
823 | { | |
824 | intersector2d.Perform(Lv, L1); | |
825 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
826 | { | |
827 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
828 | Standard_Real param1 = i2d.ParamOnFirst(); | |
829 | Standard_Real param2 = i2d.ParamOnSecond(); | |
830 | if (param1 > Precision::Confusion() && param1 < dv && | |
831 | param2 > Precision::Confusion() && param2 < d1) | |
832 | { | |
833 | continue; | |
834 | } | |
835 | } | |
836 | } | |
837 | ||
838 | // Side 2: | |
839 | if (mayIntersect(p2d2, p2d3, p2dc, p2dt)) | |
840 | { | |
841 | intersector2d.Perform(Lv, L2); | |
842 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
843 | { | |
844 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
845 | Standard_Real param1 = i2d.ParamOnFirst(); | |
846 | Standard_Real param2 = i2d.ParamOnSecond(); | |
847 | if (param1 > Precision::Confusion() && param1 < dv && | |
848 | param2 > Precision::Confusion() && param2 < d2) | |
849 | { | |
850 | continue; | |
851 | } | |
852 | } | |
853 | } | |
854 | ||
855 | // Side 3: | |
856 | if (mayIntersect(p2d3, p2d1, p2dc, p2dt)) | |
857 | { | |
858 | intersector2d.Perform(Lv, L3); | |
859 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
860 | { | |
861 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
862 | Standard_Real param1 = i2d.ParamOnFirst(); | |
863 | Standard_Real param2 = i2d.ParamOnSecond(); | |
864 | if (param1 > Precision::Confusion() && param1 < dv && | |
865 | param2 > Precision::Confusion() && param2 < d3) | |
866 | { | |
867 | continue; | |
868 | } | |
869 | } | |
870 | } | |
871 | ||
872 | //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, Standard_True); | |
873 | ||
874 | // Check intersection between the triangle & sub-voxels of the sub-voxel. | |
875 | for (Standard_Integer j = 0; j < 8; j++) | |
876 | { | |
877 | ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, j, xc, yc, zc); | |
878 | pc.SetCoord(xc, yc, zc); | |
879 | if (plane.Distance(pc) < hdiagonal4) | |
880 | { | |
881 | ElSLib::Parameters(plane, pc, uc, vc); | |
882 | p2dc.SetCoord(uc, vc); | |
883 | ||
884 | gp_Vec2d v2dct(p2dc, p2dt); | |
885 | dv = v2dct.Magnitude() - Precision::Confusion(); | |
886 | Lv.SetLocation(p2dc); | |
887 | Lv.SetDirection(v2dct); | |
888 | ||
889 | // Side 1: | |
890 | if (mayIntersect(p2d1, p2d2, p2dc, p2dt)) | |
891 | { | |
892 | intersector2d.Perform(Lv, L1); | |
893 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
894 | { | |
895 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
896 | Standard_Real param1 = i2d.ParamOnFirst(); | |
897 | Standard_Real param2 = i2d.ParamOnSecond(); | |
898 | if (param1 > Precision::Confusion() && param1 < dv && | |
899 | param2 > Precision::Confusion() && param2 < d1) | |
900 | { | |
901 | continue; | |
902 | } | |
903 | } | |
904 | } | |
905 | ||
906 | // Side 2: | |
907 | if (mayIntersect(p2d2, p2d3, p2dc, p2dt)) | |
908 | { | |
909 | intersector2d.Perform(Lv, L2); | |
910 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
911 | { | |
912 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
913 | Standard_Real param1 = i2d.ParamOnFirst(); | |
914 | Standard_Real param2 = i2d.ParamOnSecond(); | |
915 | if (param1 > Precision::Confusion() && param1 < dv && | |
916 | param2 > Precision::Confusion() && param2 < d2) | |
917 | { | |
918 | continue; | |
919 | } | |
920 | } | |
921 | } | |
922 | ||
923 | // Side 3: | |
924 | if (mayIntersect(p2d3, p2d1, p2dc, p2dt)) | |
925 | { | |
926 | intersector2d.Perform(Lv, L3); | |
927 | if (intersector2d.IsDone() && !intersector2d.ParallelElements() && intersector2d.NbPoints()) | |
928 | { | |
929 | const IntAna2d_IntPoint& i2d = intersector2d.Point(1); | |
930 | Standard_Real param1 = i2d.ParamOnFirst(); | |
931 | Standard_Real param2 = i2d.ParamOnSecond(); | |
932 | if (param1 > Precision::Confusion() && param1 < dv && | |
933 | param2 > Precision::Confusion() && param2 < d3) | |
934 | { | |
935 | continue; | |
936 | } | |
937 | } | |
938 | } | |
939 | ||
940 | ((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, j, Standard_True); | |
941 | } | |
942 | } // End of "Check level 2". | |
943 | ||
944 | } | |
945 | } // End of "Check level 1". | |
946 | ||
947 | break; | |
948 | } | |
949 | } | |
950 | } | |
951 | } | |
952 | } | |
953 | } | |
954 | } | |
f03671a0 P |
955 | |
956 | //! This macro quickly finds the min & max values among 3 variables | |
957 | #define FINDMINMAX(x0, x1, x2, min, max) \ | |
958 | min = max = x0; \ | |
959 | if(x1<min) min=x1; \ | |
960 | if(x1>max) max=x1; \ | |
961 | if(x2<min) min=x2; \ | |
962 | if(x2>max) max=x2; | |
963 | ||
964 | static bool planeBoxOverlap(const gp_Vec & normal, const double d, const gp_Pnt & maxbox) | |
965 | { | |
966 | gp_Vec vmin, vmax; | |
967 | if(normal.X() > 0.0) { vmin.SetX(-maxbox.X()); vmax.SetX(maxbox.X()); } | |
968 | else { vmin.SetX(maxbox.X()); vmax.SetX(-maxbox.X()); } | |
969 | ||
970 | if(normal.Y() > 0.0) { vmin.SetY(-maxbox.Y()); vmax.SetY(maxbox.Y()); } | |
971 | else { vmin.SetY(maxbox.Y()); vmax.SetY(-maxbox.Y()); } | |
972 | ||
973 | if(normal.Z() > 0.0) { vmin.SetZ(-maxbox.Z()); vmax.SetZ(maxbox.Z()); } | |
974 | else { vmin.SetZ(maxbox.Z()); vmax.SetZ(-maxbox.Z()); } | |
975 | ||
976 | if((normal.Dot(vmin)) + d > 0.0) return false; | |
977 | if((normal.Dot(vmax)) + d>= 0.0) return true; | |
978 | ||
979 | return false; | |
980 | } | |
981 | ||
982 | #define AXISTEST_X01(a, b, fa, fb) \ | |
983 | min = a*v0.Y() - b*v0.Z(); \ | |
984 | max = a*v2.Y() - b*v2.Z(); \ | |
985 | if(min>max) {const double tmp=max; max=min; min=tmp; } \ | |
986 | rad = fa * extents.Y() + fb * extents.Z(); \ | |
987 | if(min>rad || max<-rad) return false; | |
988 | ||
989 | #define AXISTEST_X2(a, b, fa, fb) \ | |
990 | min = a*v0.Y() - b*v0.Z(); \ | |
991 | max = a*v1.Y() - b*v1.Z(); \ | |
992 | if(min>max) {const double tmp=max; max=min; min=tmp; } \ | |
993 | rad = fa * extents.Y() + fb * extents.Z(); \ | |
994 | if(min>rad || max<-rad) return false; | |
995 | ||
996 | #define AXISTEST_Y02(a, b, fa, fb) \ | |
997 | min = b*v0.Z() - a*v0.X(); \ | |
998 | max = b*v2.Z() - a*v2.X(); \ | |
999 | if(min>max) {const double tmp=max; max=min; min=tmp; } \ | |
1000 | rad = fa * extents.X() + fb * extents.Z(); \ | |
1001 | if(min>rad || max<-rad) return false; | |
1002 | ||
1003 | #define AXISTEST_Y1(a, b, fa, fb) \ | |
1004 | min = b*v0.Z() - a*v0.X(); \ | |
1005 | max = b*v1.Z() - a*v1.X(); \ | |
1006 | if(min>max) {const double tmp=max; max=min; min=tmp; } \ | |
1007 | rad = fa * extents.X() + fb * extents.Z(); \ | |
1008 | if(min>rad || max<-rad) return false; | |
1009 | ||
1010 | #define AXISTEST_Z12(a, b, fa, fb) \ | |
1011 | min = a*v1.X() - b*v1.Y(); \ | |
1012 | max = a*v2.X() - b*v2.Y(); \ | |
1013 | if(min>max) {const double tmp=max; max=min; min=tmp; } \ | |
1014 | rad = fa * extents.X() + fb * extents.Y(); \ | |
1015 | if(min>rad || max<-rad) return false; | |
1016 | ||
1017 | #define AXISTEST_Z0(a, b, fa, fb) \ | |
1018 | min = a*v0.X() - b*v0.Y(); \ | |
1019 | max = a*v1.X() - b*v1.Y(); \ | |
1020 | if(min>max) {const double tmp=max; max=min; min=tmp; } \ | |
1021 | rad = fa * extents.X() + fb * extents.Y(); \ | |
1022 | if(min>rad || max<-rad) return false; | |
1023 | ||
1024 | // compute triangle edges | |
1025 | // - edges lazy evaluated to take advantage of early exits | |
1026 | // - fabs precomputed (half less work, possible since extents are always >0) | |
1027 | // - customized macros to take advantage of the null component | |
1028 | // - axis vector discarded, possibly saves useless movs | |
1029 | #define IMPLEMENT_CLASS3_TESTS \ | |
1030 | double rad; \ | |
1031 | double min, max; \ | |
1032 | \ | |
1033 | const double fey0 = fabs(e0.Y()); \ | |
1034 | const double fez0 = fabs(e0.Z()); \ | |
1035 | AXISTEST_X01(e0.Z(), e0.Y(), fez0, fey0); \ | |
1036 | const double fex0 = fabs(e0.X()); \ | |
1037 | AXISTEST_Y02(e0.Z(), e0.X(), fez0, fex0); \ | |
1038 | AXISTEST_Z12(e0.Y(), e0.X(), fey0, fex0); \ | |
1039 | \ | |
1040 | const double fey1 = fabs(e1.Y()); \ | |
1041 | const double fez1 = fabs(e1.Z()); \ | |
1042 | AXISTEST_X01(e1.Z(), e1.Y(), fez1, fey1); \ | |
1043 | const double fex1 = fabs(e1.X()); \ | |
1044 | AXISTEST_Y02(e1.Z(), e1.X(), fez1, fex1); \ | |
1045 | AXISTEST_Z0(e1.Y(), e1.X(), fey1, fex1); \ | |
1046 | \ | |
1047 | const gp_Vec e2 = v2 - v0; \ | |
1048 | const double fey2 = fabs(e2.Y()); \ | |
1049 | const double fez2 = fabs(e2.Z()); \ | |
1050 | AXISTEST_X2(e2.Z(), e2.Y(), fez2, fey2); \ | |
1051 | const double fex2 = fabs(e2.X()); \ | |
1052 | AXISTEST_Y1(e2.Z(), e2.X(), fez2, fex2); \ | |
1053 | AXISTEST_Z12(e2.Y(), e2.X(), fey2, fex2); | |
1054 | ||
1055 | static bool TriBoxOverlap(const gp_Pnt & p1, const gp_Pnt & p2, const gp_Pnt & p3, | |
1056 | const gp_Pnt & center, const gp_Pnt & extents) | |
1057 | { | |
1058 | // use separating axis theorem to test overlap between triangle and box | |
1059 | // need to test for overlap in these directions: | |
1060 | // 1) the {x,y,z}-directions (actually, since we use the AABB of the triangle | |
1061 | // we do not even need to test these) | |
1062 | // 2) normal of the triangle | |
1063 | // 3) crossproduct(edge from tri, {x,y,z}-directin) | |
1064 | // this gives 3x3=9 more tests | |
1065 | ||
1066 | // move everything so that the boxcenter is in (0,0,0) | |
1067 | gp_Vec v0(center, p1); | |
1068 | gp_Vec v1(center, p2); | |
1069 | gp_Vec v2(center, p3); | |
1070 | ||
1071 | // First, test overlap in the {x,y,z}-directions | |
1072 | double min,max; | |
1073 | // Find min, max of the triangle in x-direction, and test for overlap in X | |
1074 | FINDMINMAX(v0.X(), v1.X(), v2.X(), min, max); | |
1075 | if(min>extents.X() || max<-extents.X()) return false; | |
1076 | ||
1077 | // Same for Y | |
1078 | FINDMINMAX(v0.Y(), v1.Y(), v2.Y(), min, max); | |
1079 | if(min>extents.Y() || max<-extents.Y()) return false; | |
1080 | ||
1081 | // Same for Z | |
1082 | FINDMINMAX(v0.Z(), v1.Z(), v2.Z(), min, max); | |
1083 | if(min>extents.Z() || max<-extents.Z()) return false; | |
1084 | ||
1085 | // 2) Test if the box intersects the plane of the triangle | |
1086 | // compute plane equation of triangle: normal*x+d=0 | |
1087 | // ### could be precomputed since we use the same leaf triangle several times | |
1088 | const gp_Vec e0 = v1 - v0; | |
1089 | const gp_Vec e1 = v2 - v1; | |
1090 | const gp_Vec normal = e0.Crossed(e1); | |
1091 | const double d = -normal.Dot(v0); | |
1092 | if(!planeBoxOverlap(normal, d, extents)) return false; | |
1093 | ||
1094 | // 3) "Class III" tests | |
1095 | //if(mFullPrimBoxTest) | |
1096 | { | |
1097 | IMPLEMENT_CLASS3_TESTS | |
1098 | } | |
1099 | ||
1100 | return true; | |
1101 | } | |
1102 | ||
1103 | void Voxel_FastConverter::ComputeVoxelsNearTriangle(const gp_Pnt& p1, | |
1104 | const gp_Pnt& p2, | |
1105 | const gp_Pnt& p3, | |
1106 | const gp_Pnt& extents, | |
1107 | const gp_Pnt& extents2, | |
1108 | const gp_Pnt& extents4, | |
1109 | const Standard_Integer ixmin, | |
1110 | const Standard_Integer iymin, | |
1111 | const Standard_Integer izmin, | |
1112 | const Standard_Integer ixmax, | |
1113 | const Standard_Integer iymax, | |
1114 | const Standard_Integer izmax) const | |
1115 | { | |
1116 | gp_Pnt pc; | |
1117 | Standard_Real xc, yc, zc; | |
1118 | Standard_Integer ix, iy, iz; | |
1119 | ||
1120 | Voxel_DS* ds = (Voxel_DS*) myVoxels; | |
1121 | for (ix = ixmin; ix <= ixmax; ix++) | |
1122 | { | |
1123 | for (iy = iymin; iy <= iymax; iy++) | |
1124 | { | |
1125 | for (iz = izmin; iz <= izmax; iz++) | |
1126 | { | |
1127 | ds->GetCenter(ix, iy, iz, xc, yc, zc); | |
1128 | pc.SetCoord(xc, yc, zc); | |
1129 | ||
1130 | if(TriBoxOverlap(p1, p2, p3, pc, extents)) | |
1131 | { | |
1132 | // Set positive value to this voxel: | |
1133 | switch (myIsBool) | |
1134 | { | |
1135 | case 0: | |
1136 | ((Voxel_ColorDS*) myVoxels)->Set(ix, iy, iz, 15); | |
1137 | break; | |
1138 | case 1: | |
1139 | ((Voxel_BoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True); | |
1140 | break; | |
1141 | case 2: | |
1142 | { | |
1143 | //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, Standard_True); | |
1144 | ||
1145 | // Check intersection between the triangle & sub-voxels of the voxel. | |
1146 | for (Standard_Integer i = 0; i < 8; i++) | |
1147 | { | |
1148 | ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, xc, yc, zc); | |
1149 | pc.SetCoord(xc, yc, zc); | |
1150 | if(TriBoxOverlap(p1, p2, p3, pc, extents2)) | |
1151 | { | |
1152 | //((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, Standard_True); | |
1153 | ||
1154 | // Check intersection between the triangle & sub-voxels of the sub-voxel. | |
1155 | for (Standard_Integer j = 0; j < 8; j++) | |
1156 | { | |
1157 | ((Voxel_ROctBoolDS*) myVoxels)->GetCenter(ix, iy, iz, i, j, xc, yc, zc); | |
1158 | pc.SetCoord(xc, yc, zc); | |
1159 | if(TriBoxOverlap(p1, p2, p3, pc, extents4)) | |
1160 | { | |
1161 | ((Voxel_ROctBoolDS*) myVoxels)->Set(ix, iy, iz, i, j, Standard_True); | |
1162 | } | |
1163 | } // End of "Check level 2". | |
1164 | ||
1165 | } | |
1166 | } // End of "Check level 1". | |
1167 | break; | |
1168 | } | |
1169 | } | |
1170 | } | |
1171 | } | |
1172 | } | |
1173 | } | |
1174 | } |