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