projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024057: Eliminate compiler warning C4100 in MSVC++ with warning level 4
[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>
4bee43a9 38#include <BRepClass3d_SolidClassifier.hxx>
7fd59977 39
40// Printing the progress in stdout.
41//#define CONV_DUMP
42
43Voxel_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),
54 myNbThreads(nbthreads),myIsBool(2),
c5e9fb8b
P		 55 myNbTriangles(0),
56 myUseExistingTriangulation(useExistingTriangulation)
7fd59977 57{
58 Init();
59}
60
61Voxel_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),
72 myNbThreads(nbthreads),myIsBool(1),
c5e9fb8b
P		 73 myNbTriangles(0),
74 myUseExistingTriangulation(useExistingTriangulation)
7fd59977 75{
76 Init();
77}
78
79Voxel_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),
90 myNbThreads(nbthreads),myIsBool(0),
c5e9fb8b
P		 91 myNbTriangles(0),
92 myUseExistingTriangulation(useExistingTriangulation)
7fd59977 93{
94 Init();
95}
96
97void 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
162void Voxel_FastConverter::Destroy()
163{
164
165}
166
167Standard_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		 304Standard_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 439Standard_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 535Standard_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 632void 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
678static 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
692void 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
964static 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
1055static 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
1103void 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}
Open CASCADE Technology repository