Производительность булевых OpenCascade

В недавно вышедшем номере журнала «САПР и Графика» обнаружилась любопытная информация о том, что булевы операции OpenCascade работают довольно медленно. Достаточно медленно, чтобы отказаться от этой библиотеки. Вот картинка, которая выглядит как приговор нашему старому доброму каскейду:

Однако есть нюанс. Мы видим, что тестировалась версия 6.5.3, которая была опубликована 24 апреля 2012 года. Считай, добрых пять лет назад. Можно кивать на то, что область геометрического моделирования, дескать, консервативна, да и что могло поменяться в библиотеке, которая, как известно, разрабатывается в «вечно замороженном углу глобуса». Но с другой стороны, как мы опять-таки знаем, наибольшее количество багфиксов OCCT относится именно к алгоритмам геометрического моделирования. И за пять лет произошло очень много любопытных вещей. Например, были радикально улучшены булевы операции. Пренебрежительное отношение к нашей библиотеке, с которым подчас приходится сталкиваться, нередко проистекает из стереотипа мышления, дескать, «бесплатно значит плохо». Не будем спорить. Будем показывать как на самом деле.

Оговорюсь сразу, что никаких сранений каскейда с конкурирующими ядрами не будет. Такое сравнение возможно только при наличии хорошего опыта работы со всеми тестируемыми библиотеками. У меня такого опыта нет и, откровенно говоря, мало у кого он есть. Поэтому многие «объективные» сравнения при внимательном рассмотрении оказываются умышленной или неосознанной спекуляцией. Сначала у меня возникло желание сравнить две версии OpenCascade, а именно, опозорившуюся 6.5.3 и свежий (уже почти выпущенный официально) релиз 7.2. Однако и это желание ваш покорный слуга подавил. Зачем валять дурака? Пусть читатель посмотрит на результаты актуального эксперимента и ответит сам, устраивает его увиденное или нет.

Строим шестиугольный примитив

Для построения шестиугольного контура воспользуемся точной формулой, которую можно найти, например, здесь. Код с использованием OpenCascade выглядит примерно следующим образом:

// Get hexagon poles in 2D
gp_XY Poles2d[6];
HexagonPoles(gp::Origin2d().XY(),
             1.0,
             Poles2d[0],
             Poles2d[1],
             Poles2d[2],
             Poles2d[3],
             Poles2d[4],
             Poles2d[5]);
  
// Choose working plane
Handle(Geom_Plane) workPlane = new Geom_Plane( gp_Pln( gp::Origin(), gp::DZ() ) );
  
// Populate wire builder
BRepBuilderAPI_MakePolygon mkPoly;
for ( int k = 0; k < 6; ++k )
  mkPoly.Add( workPlane->Value( Poles2d[k].X(), Poles2d[k].Y() ) );
//
mkPoly.Add( workPlane->Value( Poles2d[0].X(), Poles2d[0].Y() ) );
  
// Build wire
mkPoly.Build();
const TopoDS_Wire& wire = mkPoly.Wire();

Функция HexagonPoles реализована так:

void HexagonPoles(const gp_XY& center,
                  const double dist2Pole,
                  gp_XY&       P1,
                  gp_XY&       P2,
                  gp_XY&       P3,
                  gp_XY&       P4,
                  gp_XY&       P5,
                  gp_XY&       P6)
{
  std::vector<gp_XY*> res = {&P1, &P2, &P3, &P4, &P5, &P6};
  
  // Taken from http://www.redblobgames.com/grids/hexagons/
  for ( int i = 0; i < 6; ++i )
  {
    const double angle_deg = 60*i + 30;
    const double angle_rad = angle_deg / 180.0 * M_PI;
    
    const double x = center.X() + dist2Pole * Cos(angle_rad);
    const double y = center.X() + dist2Pole * Sin(angle_rad);
    
    res[i]->SetX(x);
    res[i]->SetY(y);
  }
}

Построение грани и призмы выполняется тривиально:

// Build face
TopoDS_Face face = BRepBuilderAPI_MakeFace(wire);
  
// Build prism
TopoDS_Shape
  prism = BRepPrimAPI_MakePrism(face,
                                gp_Vec( workPlane->Axis().Direction() )*0.25);

Строим цилиндры

Подозреваю, что цилиндры из статьи моделируют что-то вроде стержней в активной зоне реактора и, конечно, не могут располагаться как попало. Мы, однако, с технологической подоплекой не знакомы. Поэтому будем создавать примитивы так, чтобы они просто-напросто попали в шестиугольную пластинку. Начать можно с того, чтобы создать избыточное количество цилиндров в объеме, ограничивающем пластину.

Чтобы отфильтровать цилиндры, выпадающие из материала пластины, мы можем использовать классификатор IntTools_FClass2d для граней. Центр каждого цилиндра проверяем на принадлежность базовой грани. Если он выпадает (классификатор возвращает статус OUT), то такой цилиндр выбраковывается.

// Build bounding box for the hexa
double xMin, yMin, zMin, xMax, yMax, zMax;
asiAlgo_Utils::Bounds(prism, xMin, yMin, zMin, xMax, yMax, zMax);
  
const double r = 0.02;
const double h = 0.25;
const double step = 0.06;
  
BRep_Builder BB;
TopoDS_Compound cyls;
BB.MakeCompound(cyls);
  
// List for cutter
TopTools_ListOfShape tools;
  
// Prepare two-dimensional classifier
IntTools_FClass2d classifier;
classifier.Init( face, Precision::Confusion() );
  
// Count number of primitives which pass the test
int cylCount = 0;
  
double xcurr = xMin;
do
{
  double ycurr = yMin;
  
  do
  {
    // Classify center point
    const TopAbs_State state = classifier.Perform( gp_Pnt2d(xcurr, ycurr) );
    //
    if ( state == TopAbs_IN )
    {
      gp_Pnt cylcenter(xcurr, ycurr, 0.0);
      TopoDS_Shape cyl = BRepPrimAPI_MakeCylinder( gp_Ax2(cylcenter, gp::DZ() ), r, 0.25 );
      //
      BB.Add(cyls, cyl);
      tools.Append(cyl);
      
      cylCount++;
    }
    
    ycurr += step;
  }
  while ( ycurr < yMax );
  
  xcurr += step;
}
while ( xcurr < xMax );

К сожалению, простая классификация не дает нужного эффекта. Центры некоторых цилиндров формально принадлежат внутренности базовой грани, но нас они не устраивают, поскольку выпадают на границу.

Чтобы улучшить результат, вместо метода Perform() классификатора IntTools_FClass2d будем использовать другой метод, принимающий в качестве дополнительного агрумента ДОПУСК:

// Classify center point
const TopAbs_State state = classifier.TestOnRestriction( gp_Pnt2d(xcurr, ycurr), 0.05 );

Булево вычитание

Для булева вычитания используем следующий код:

//! Performs Boolean Cut of a tool shape from the given object shape.
//! param Object [in] object to cut from.
//! param Tool   [in] tool to cut out.
//! param fuzzy  [in] fuzzy tolerance.
//! eturn result shape.
TopoDS_Shape BooleanCut(const TopoDS_Shape& Object,
                        const TopoDS_Shape& Tool,
                        const double        fuzzy)
{
  // Prepare the arguments
  TopTools_ListOfShape BOP_args;
  BOP_args.Append(Object);
  BOP_args.Append(Tool);
  
  // Prepare data structure (calculate interferences)
  Handle(NCollection_BaseAllocator) Alloc = new NCollection_IncAllocator;
  //
  BOPAlgo_PaveFiller DSFiller(Alloc);
  DSFiller.SetArguments(BOP_args);
  DSFiller.SetRunParallel(0);
  DSFiller.SetFuzzyValue(fuzzy);
  DSFiller.Perform();
  
  // Check data structure
  bool hasErr = DSFiller.HasErrors();
  if ( hasErr )
  {
    std::cout << "Error: cannot cut" << std::endl;
    return TopoDS_Shape();
  }
  
  // Run BOP
  BOPAlgo_BOP BOP(Alloc);
  BOP.AddArgument(Object);
  BOP.AddTool(Tool);
  BOP.SetRunParallel(0);
  BOP.SetOperation(BOPAlgo_CUT);
  BOP.PerformWithFiller(DSFiller);
  hasErr = BOP.HasErrors();
  if ( hasErr )
  {
    std::cout << "Error: cannot cut the part model from the stock model" << std::endl;
    return TopoDS_Shape();
  }
  
  return BOP.Shape();
}

Итоги

Может быть у меня слишком мощное железо?

Обновление от 2018-03-24

Тов. Михаил (см. комментарии ниже) справедливо заметил, что вместо вычитания трехмерных тел, для получения нужного результата достаточно вычесть кружки из шестиугольного эскиза, после чего сделать операцию вытягивания. Результат такого эксперимента порадовал: построение тела ускоряется примерно в 10 раз. В связи с этим привожу итоговую таблицу: