zippy84 / vtkbool

A new boolean operations filter for VTK
Apache License 2.0
166 stars 40 forks source link

Boolean union of touching objects returns empty mesh #80

Open mauigna06 opened 4 months ago

mauigna06 commented 4 months ago

Hi,

Thank you for vtkbool, it is great :)

I have a bug report

this is my starting model: image

then I add this one: image

Till there behavior is correct.

But then when trying to merge this one too: image

And then the result is an empty mesh


Here are the objects in vtk format: object0.zip object1.zip object2.zip


By the way, the operation is successful if I apply some very small translation (e.g. 0.1mm) to some of the objects in some random direction

mauigna06 commented 4 months ago

Also, another question? What's better? to do 3 different union operations, or first do the union of object1 and object2 (that have less points) and then join the result to object0? (robustness-wise and minimizing-execution-time wise)

zippy84 commented 4 months ago

@mauigna06 Thanks for your report. I tried to combine it with Python:

#!/usr/bin/env python
# *-* coding: UTF-8 *-*

import sys
sys.path.extend(['/home/zippy/vtkbool/build/lib/python3.12/site-packages/vtkbool'])

from vtkmodules.vtkIOLegacy import vtkPolyDataReader, vtkPolyDataWriter
from vtkmodules.vtkIOGeometry import vtkSTLReader, vtkSTLWriter
from vtkmodules.vtkFiltersCore import vtkCleanPolyData, vtkPolyDataNormals, vtkTriangleFilter
from vtkBool import vtkPolyDataBooleanFilter

rA = vtkPolyDataReader()
rA.SetFileName('object0.vtk')

rB = vtkPolyDataReader()
rB.SetFileName('object1.vtk')

rC = vtkPolyDataReader()
rC.SetFileName('object2.vtk')

bfA = vtkPolyDataBooleanFilter()
bfA.SetInputConnection(0, rA.GetOutputPort())
bfA.SetInputConnection(1, rB.GetOutputPort())

bfB = vtkPolyDataBooleanFilter()
bfB.SetInputConnection(0, bfA.GetOutputPort())
bfB.SetInputConnection(1, rC.GetOutputPort())

writerA = vtkPolyDataWriter()
writerA.SetFileName('unionA.vtk')
writerA.SetInputConnection(bfA.GetOutputPort())

clean = vtkCleanPolyData()
clean.SetInputConnection(bfB.GetOutputPort())

writerB = vtkPolyDataWriter()
writerB.SetFileName('unionB.vtk')
writerB.SetInputConnection(clean.GetOutputPort())
writerB.Update()

The result is valid. I got no errors.

Here is what I get: unionB.zip

I think you are not using the current version. 3.0.1 has a few important bugfixes.

mauigna06 commented 4 months ago

Hi @zippy84. Thanks a lot for your answer. I'll share the models to reproduce the problem on latest vtkbool version

mauigna06 commented 4 months ago

Hi @zippy84. Thank for your great work on this tool.

Download test data for current code - 85aac4470d6740c0beb7d48b3f329ffb5df188cf - bugs

Context:

Reproduce bugs:

Thanks a lot Ronald :)

zippy84 commented 3 months ago

@mauigna06 Thanks for sharing your data. I figured out, why meshC - meshD leads to that error. There is a missing line between polys 142830 and 47968 of newPdA and newPdB. There is a rounding issue here: https://github.com/zippy84/vtkbool/blob/85aac4470d6740c0beb7d48b3f329ffb5df188cf/vtkPolyDataContactFilter.cxx#L506-L507

Here is what happens:

#include <string>

int main() {
  // InterPolys() -> idA 142830, idB 47968 
  // pt 2012
  double a = -110.6449891, b = -110.6449823;

  // wrong result
  std::cout << std::lround(a*1e5) << ", " << std::lround(b*1e5) << '\n'; // -> -11064499, -11064498

  // correct
  std::cout << static_cast<long>(a*1e5) << ", " << static_cast<long>(b*1e5) << '\n'; // -11064498, -11064498
}

If you change std::lround to _staticcast, you will get the correct difference:

grafik

I don't know if this is a bug and if I should fix it?!

mauigna06 commented 3 months ago

Hi @zippy84

I don't know if this is a bug and if I should fix it?!

I don't know either, I think there is no conclusive answer:

https://www.perplexity.ai/search/is-it-better-to-compare-int-th-dhdZFkf.RcS41BNRmo.u3g#1

I compiled the code below and it works


I asked Copilot, the problem is on the lines copied below, solve it:

const long a = std::lround(lhs.t * 1e5),
    b = std::lround(rhs.t * 1e5);

return a < b;

It said:

The problem with the lines in question is related to the precision and potential rounding errors when comparing floating-point numbers after scaling and rounding them. To solve this, we can improve the comparison by considering a tolerance level for the floating-point comparison, which can help in addressing the issue of comparing very close values that should be considered equal. This approach avoids the pitfalls of direct floating-point comparisons and the issues introduced by rounding.

Here's how you can solve it:

  1. Define a tolerance level for comparison. This value represents the smallest difference between lhs.t and rhs.t that we consider significant.
  2. Compare the floating-point numbers lhs.t and rhs.t directly, considering the tolerance to decide if they are effectively equal or if one is less than the other.

Here's the revised code:

const double tolerance = 1e-5; // Define a tolerance level for comparison

// Calculate the difference once to avoid doing it multiple times
const double diff = lhs.t - rhs.t;

// Check if the difference is within the tolerance level
if (std::fabs(diff) < tolerance) {
    return false; // Consider them equal, so lhs is not less than rhs
} else {
    return diff < 0; // If diff is negative, lhs is less than rhs
}

This solution directly compares the t values of lhs and rhs, using a tolerance to handle the precision issues inherent in floating-point arithmetic. This method is more accurate and reliable for floating-point comparisons than scaling and rounding, especially when dealing with very small differences between numbers.


Now, meshC - meshD succeeds (in 2.7 seconds)