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