GCC Code Coverage Report

Directory:	./
File:	delaunay/LFS.cpp
Date:	2026-09-27 03:24:14

	Exec	Total	Coverage
Lines:	139	140	99.3%
Functions:	5	5	100.0%
Branches:	171	286	59.8%

  
      Line
      Branch
      Exec
      Source
    
      /*
    
       *  Copyright (c) 2000-2022 Inria
    
       *  All rights reserved.
    
       *
    
       *  Redistribution and use in source and binary forms, with or without
    
       *  modification, are permitted provided that the following conditions are met:
    
       *
    
       *  * Redistributions of source code must retain the above copyright notice,
    
       *  this list of conditions and the following disclaimer.
    
       *  * Redistributions in binary form must reproduce the above copyright notice,
    
       *  this list of conditions and the following disclaimer in the documentation
    
       *  and/or other materials provided with the distribution.
    
       *  * Neither the name of the ALICE Project-Team nor the names of its
    
       *  contributors may be used to endorse or promote products derived from this
    
       *  software without specific prior written permission.
    
       *
    
       *  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
    
       *  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
    
       *  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
    
       *  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
    
       *  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
    
       *  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
    
       *  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
    
       *  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
    
       *  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
    
       *  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
    
       *  POSSIBILITY OF SUCH DAMAGE.
    
       *
    
       *  Contact: Bruno Levy
    
       *
    
       *     https://www.inria.fr/fr/bruno-levy
    
       *
    
       *     Inria,
    
       *     Domaine de Voluceau,
    
       *     78150 Le Chesnay - Rocquencourt
    
       *     FRANCE
    
       *
    
       */
    
      #include <geogram/delaunay/LFS.h>
    
      #include <geogram/basic/geometry.h>
    
      #include <geogram/basic/logger.h>
    
      #include <geogram/basic/command_line.h>
    
      namespace {
    
          using namespace GEO;
    
          const index_t pid_size = 6;
    
          int pid[pid_size][4] = {
    
              {0, 1, 2, 3},
    
              {0, 2, 1, 3},
    
              {0, 3, 1, 2},
    
              {1, 2, 0, 3},
    
              {1, 3, 2, 0},
    
              {2, 3, 0, 1}
    
          };
    
          /**
    
           * \brief Computes the circumcenter and squared radius of a tetrahedron.
    
           * \param[in] p first vertex of the tetrahedron
    
           * \param[in] q second vertex of the tetrahedron
    
           * \param[in] r third vertex of the tetrahedron
    
           * \param[in] s fourth vertex of the tetrahedron
    
           * \param[out] t_circumcenter computed circumcenter of the tetrahedron
    
           * \param[out] squared_radius computed squared radius of the tetrahedron
    
           * \param[out] dihedral_angles if non-nullptr,
    
           *   computed dihedral angles of the tetrahedron
    
           * \return true if successful, false otherwise
    
           *  (for instance, if the tetrahedron is flat).
    
           */
    
      55451
          bool tetra_circumcenter_squaredradius(
    
              const vec3& p, const vec3& q,
    
              const vec3& r, const vec3& s,
    
              vec3& t_circumcenter, double& squared_radius,
    
              double* dihedral_angles = nullptr
    
          ) {
    
      55451
              vec3 qp = q - p;
    
      55451
              vec3 rp = r - p;
    
      55451
              vec3 sp = s - p;
    
      55451
              double tet_vol_x_12 = 2.0 * dot(qp, cross(rp, sp));
    
              // it is not a safe check, one should avoid
    
              // to input a degenerated tetra.
    
        1/2✗ Branch 0 not taken.
✓ Branch 1 taken 55451 times.

      55451
              if(tet_vol_x_12 == 0.0) {
    
      ✗
                  return false;
    
              }
    
              vec3 d_t_circumcenter =
    
        1/2✓ Branch 2 taken 55451 times.
✗ Branch 3 not taken.

      55451
                  length2(qp) * cross(rp, sp) +
    
        1/2✓ Branch 2 taken 55451 times.
✗ Branch 3 not taken.

      110902
                  length2(rp) * cross(sp, qp) +
    
        1/2✓ Branch 2 taken 55451 times.
✗ Branch 3 not taken.

      110902
                  length2(sp) * cross(qp, rp);
    
        1/2✓ Branch 1 taken 55451 times.
✗ Branch 2 not taken.

      55451
              squared_radius = length2(
    
      55451
                  d_t_circumcenter) / (tet_vol_x_12 * tet_vol_x_12
    
                                      );
    
      55451
              d_t_circumcenter = (1.0 / tet_vol_x_12) * d_t_circumcenter;
    
      55451
              t_circumcenter = p + d_t_circumcenter;
    
        1/2✓ Branch 0 taken 55451 times.
✗ Branch 1 not taken.

      55451
              if(dihedral_angles) {
    
                  // compute Dihedral angle directly
    
                  // the following computation is not optimized.
    
        2/2✓ Branch 1 taken 221804 times.
✓ Branch 2 taken 55451 times.

      277255
                  vec3 point[4];
    
      55451
                  point[0] = p;
    
      55451
                  point[1] = q;
    
      55451
                  point[2] = r;
    
      55451
                  point[3] = s;
    
        2/2✓ Branch 0 taken 332706 times.
✓ Branch 1 taken 55451 times.

      388157
                  for(index_t i = 0; i < pid_size; i++) {
    
      332706
                      vec3 b1 = point[pid[i][0]] - point[pid[i][2]];
    
      332706
                      vec3 b2 = point[pid[i][1]] - point[pid[i][0]];
    
      332706
                      vec3 b3 = point[pid[i][3]] - point[pid[i][1]];
    
      332706
                      vec3 b2b3 = cross(b2, b3);
    
      332706
                      dihedral_angles[i] = ::fabs(
    
      332706
                          ::atan2(
    
        1/2✓ Branch 1 taken 332706 times.
✗ Branch 2 not taken.

      332706
                              length(b2) * dot(b1, b2b3),
    
      665412
                              dot(cross(b1, b2), b2b3)
    
                          )
    
                      );
    
                  }
    
              }
    
      55451
              return true;
    
          }
    
          /**
    
           * \brief Gets a Delaunay vertex by global index.
    
           * \param[in] delaunay the Delaunay triangulation
    
           * \param[in] v the index of the vertex
    
           * \return a const reference to the vertex, as a vec3
    
           */
    
      633474
          inline const vec3& delaunay_vertex(Delaunay* delaunay, index_t v) {
    
      633474
              return *(const vec3*) delaunay->vertex_ptr(v);
    
          }
    
          /**
    
           * \brief Gets a Delaunay vertex by tetrahedron index
    
           *  and local vertex index.
    
           * \param[in] delaunay the Delaunay triangulation
    
           * \param[in] c the index of the tetrahedron
    
           * \param[in] lv the local index of the vertex (0,1,2 or 3)
    
           *  in tetrahedron \p c.
    
           * \return a const reference to the vertex, as a vec3.
    
           */
    
      221804
          inline const vec3& delaunay_tet_vertex(
    
              Delaunay* delaunay, index_t c, index_t lv
    
          ) {
    
      221804
              return delaunay_vertex(delaunay, index_t(delaunay->cell_vertex(c, lv)));
    
          }
    
          index_t tet_facet_vertex[4][3] = {
    
              {1, 2, 3},
    
              {0, 3, 2},
    
              {3, 0, 1},
    
              {2, 1, 0}
    
          };
    
          /**
    
           * \brief Computes the normal to a tetrahedron facet.
    
           * \param[in] delaunay the Delaunay triangulation
    
           * \param[in] t the index of the tetrahedron
    
           * \param[in] f the local index (0,1,2 or 3) of
    
           *  the facet in the tetrahedron \p t
    
           * \return the normal to the facet \p f or tetrahedron \p t
    
           */
    
      6250
          vec3 delaunay_facet_normal(
    
              Delaunay* delaunay, index_t t, index_t f
    
          ) {
    
        1/6✗ Branch 0 not taken.
✓ Branch 1 taken 6250 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.

      6250
              geo_debug_assert(f < 4);
    
      6250
              index_t v1 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][0]));
    
      6250
              index_t v2 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][1]));
    
      6250
              index_t v3 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][2]));
    
      6250
              const vec3& p1 = delaunay_vertex(delaunay, v1);
    
      6250
              const vec3& p2 = delaunay_vertex(delaunay, v2);
    
      6250
              const vec3& p3 = delaunay_vertex(delaunay, v3);
    
      6250
              return cross(p2 - p1, p3 - p1);
    
          }
    
      }
    
      /****************************************************************************/
    
      namespace GEO {
    
      1
          void LocalFeatureSize::init(index_t nb_pts, const double* pts) {
    
              // Note: I need PDEL here instead of ANN/BNN since I need
    
              //  to access the cells (and ANN/BNN only give me the
    
              //  neighbors)
    
        4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.

      2
              Logger::out("LFS") << "Delaunay" << std::endl;
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.

      1
              Delaunay_var delaunay = Delaunay::create(3, "PDEL");
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.

      1
              delaunay->set_vertices(nb_pts, pts);
    
        4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.

      1
              Logger::out("LFS") << "Done Delaunay" << std::endl;
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.

      1
              vector<vec3> circumcenter(delaunay->nb_cells());
    
              vector<bool> voronoi_cell_is_infinite(
    
                  delaunay->nb_vertices(), false
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.

      1
              );
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.

      1
              vector<bool> is_sliver(delaunay->nb_cells());
    
              static constexpr index_t NO_POLE = NO_INDEX;
    
              // Index of incident tet whose circumcenter is the pole
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.

      1
              vector<index_t> positive_pole(delaunay->nb_vertices(),NO_POLE);
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.

      1
              vector<index_t> negative_pole(delaunay->nb_vertices(),NO_POLE);
    
              vector<vec3> avg_infinite_dir(
    
      1
                  delaunay->nb_vertices(), vec3(0, 0, 0)
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.

      1
              );
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.

      1
              vector<double> dist(delaunay->nb_vertices(), 0.0);
    
        4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.

      1
              Logger::out("LFS") << "(1) Circumcenters and slivers" << std::endl;
    
              // Step 1: compute circumcenters and check for slivers
    
      1
              const double sliver_quality = sliver_angle_threshold_ / 180.0 * M_PI;
    
        3/4✓ Branch 1 taken 55452 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55451 times.
✓ Branch 5 taken 1 times.

      55452
              for(index_t t = 0; t < delaunay->nb_cells(); t++) {
    
        1/2✓ Branch 2 taken 55451 times.
✗ Branch 3 not taken.

      55451
                  const vec3& p = delaunay_tet_vertex(delaunay, t, 0);
    
        1/2✓ Branch 2 taken 55451 times.
✗ Branch 3 not taken.

      55451
                  const vec3& q = delaunay_tet_vertex(delaunay, t, 1);
    
        1/2✓ Branch 2 taken 55451 times.
✗ Branch 3 not taken.

      55451
                  const vec3& r = delaunay_tet_vertex(delaunay, t, 2);
    
        1/2✓ Branch 2 taken 55451 times.
✗ Branch 3 not taken.

      55451
                  const vec3& s = delaunay_tet_vertex(delaunay, t, 3);
    
                  double dihedral_angle[6];
    
                  double tet_radius;
    
      110902
                  bool ok = tetra_circumcenter_squaredradius(
    
                      p, q, r, s,
    
        2/4✓ Branch 1 taken 55451 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55451 times.
✗ Branch 5 not taken.

      55451
                      circumcenter[t], tet_radius, dihedral_angle
    
                  );
    
        2/2✓ Branch 0 taken 332449 times.
✓ Branch 1 taken 55375 times.

      387824
                  for(index_t a = 0; a < 6; a++) {
    
        2/2✓ Branch 0 taken 76 times.
✓ Branch 1 taken 332373 times.

      332449
                      if(!ok) {
    
      76
                          break;
    
                      }
    
      332373
                      ok = ok &&
    
        3/4✓ Branch 0 taken 332373 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 332295 times.
✓ Branch 3 taken 78 times.

      664668
                          dihedral_angle[a] >= sliver_quality &&
    
        2/2✓ Branch 0 taken 332287 times.
✓ Branch 1 taken 8 times.

      332295
                          dihedral_angle[a] <= M_PI - sliver_quality;
    
                  }
    
      55451
                  is_sliver[t] = !ok;
    
              }
    
        4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.

      1
              Logger::out("LFS") << "(2) Positive poles" << std::endl;
    
              // Step 2: compute positive poles
    
        3/4✓ Branch 1 taken 55452 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55451 times.
✓ Branch 5 taken 1 times.

      55452
              for(index_t t = 0; t < delaunay->nb_cells(); t++) {
    
        2/2✓ Branch 2 taken 86 times.
✓ Branch 3 taken 55365 times.

      55451
                  if(is_sliver[t]) {
    
      86
                      continue;
    
                  }
    
      55365
                  index_t f_inf = NO_INDEX;
    
        2/2✓ Branch 0 taken 209269 times.
✓ Branch 1 taken 49115 times.

      258384
                  for(index_t f = 0; f < 4; f++) {
    
        4/6✓ Branch 1 taken 209269 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 209269 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6250 times.
✓ Branch 7 taken 203019 times.

      209269
                      if(delaunay->cell_adjacent(t, f) == NO_INDEX) {
    
      6250
                          f_inf = index_t(f);
    
      6250
                          break;
    
                      }
    
                  }
    
        2/2✓ Branch 0 taken 49115 times.
✓ Branch 1 taken 6250 times.

      55365
                  if(f_inf == NO_INDEX) {
    
                      // tet t does not have facet on border
    
        1/2✓ Branch 1 taken 49115 times.
✗ Branch 2 not taken.

      49115
                      const vec3& c = circumcenter[t];
    
        2/2✓ Branch 0 taken 196460 times.
✓ Branch 1 taken 49115 times.

      245575
                      for(index_t lv = 0; lv < 4; lv++) {
    
        2/4✓ Branch 1 taken 196460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 196460 times.
✗ Branch 5 not taken.

      196460
                          index_t iv = index_t(delaunay->cell_vertex(t, lv));
    
        1/2✓ Branch 2 taken 196460 times.
✗ Branch 3 not taken.

      196460
                          const vec3& v = delaunay_vertex(delaunay, iv);
    
        1/2✓ Branch 2 taken 196460 times.
✗ Branch 3 not taken.

      196460
                          double d = length2(c - v);
    
        3/4✓ Branch 1 taken 196460 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 34779 times.
✓ Branch 4 taken 161681 times.

      196460
                          if(d > dist[iv]) {
    
        1/2✓ Branch 1 taken 34779 times.
✗ Branch 2 not taken.

      34779
                              dist[iv] = d;
    
        1/2✓ Branch 1 taken 34779 times.
✗ Branch 2 not taken.

      34779
                              positive_pole[iv] = t;
    
                          }
    
                      }
    
                  } else {
    
                      // tet t has facet f_inf on border
    
                      // positive pole = average direction of infinite Voronoi edges
    
        2/4✓ Branch 1 taken 6250 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6250 times.
✗ Branch 5 not taken.

      6250
                      index_t v = index_t(delaunay->cell_vertex(t, index_t(f_inf)));
    
      6250
                      voronoi_cell_is_infinite[v] = true;
    
        1/2✓ Branch 1 taken 6250 times.
✗ Branch 2 not taken.

      6250
                      avg_infinite_dir[v] +=
    
        1/2✓ Branch 2 taken 6250 times.
✗ Branch 3 not taken.

      12500
                          delaunay_facet_normal(delaunay, t, index_t(f_inf));
    
                  }
    
              }
    
        4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.

      1
              Logger::out("LFS") << "(3) Negative poles" << std::endl;
    
              // Step 3: compute negative poles
    
      1
              std::fill(dist.begin(), dist.end(), 0.0);
    
        3/4✓ Branch 1 taken 55452 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55451 times.
✓ Branch 5 taken 1 times.

      55452
              for(index_t t = 0; t < delaunay->nb_cells(); t++) {
    
        2/2✓ Branch 2 taken 86 times.
✓ Branch 3 taken 55365 times.

      55451
                  if(is_sliver[t]) {
    
      86
                      continue;
    
                  }
    
      55365
                  bool t_is_infinite = false;
    
        2/2✓ Branch 0 taken 209269 times.
✓ Branch 1 taken 49115 times.

      258384
                  for(index_t f = 0; f < 4; f++) {
    
        4/6✓ Branch 1 taken 209269 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 209269 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 6250 times.
✓ Branch 7 taken 203019 times.

      209269
                      if(delaunay->cell_adjacent(t, f) == NO_INDEX) {
    
      6250
                          t_is_infinite = true;
    
      6250
                          break;
    
                      }
    
                  }
    
        2/2✓ Branch 0 taken 6250 times.
✓ Branch 1 taken 49115 times.

      55365
                  if(t_is_infinite) {
    
      6250
                      continue;
    
                  }
    
        1/2✓ Branch 1 taken 49115 times.
✗ Branch 2 not taken.

      49115
                  const vec3& c = circumcenter[t];
    
        2/2✓ Branch 0 taken 196460 times.
✓ Branch 1 taken 49115 times.

      245575
                  for(index_t lv = 0; lv < 4; lv++) {
    
        2/4✓ Branch 1 taken 196460 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 196460 times.
✗ Branch 5 not taken.

      196460
                      index_t iv = index_t(delaunay->cell_vertex(t, lv));
    
        1/2✓ Branch 2 taken 196460 times.
✗ Branch 3 not taken.

      196460
                      const vec3& v = delaunay_vertex(delaunay, iv);
    
      196460
                      vec3 N;
    
        2/2✓ Branch 2 taken 40411 times.
✓ Branch 3 taken 156049 times.

      196460
                      if(voronoi_cell_is_infinite[iv]) {
    
        1/2✓ Branch 1 taken 40411 times.
✗ Branch 2 not taken.

      40411
                          N = -avg_infinite_dir[iv];
    
                      } else {
    
        2/4✓ Branch 1 taken 156049 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 156049 times.
✗ Branch 5 not taken.

      156049
                          N = v - circumcenter[positive_pole[iv]];
    
                      }
    
      196460
                      double d = dot(c - v, N);
    
        3/4✓ Branch 1 taken 196460 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 25805 times.
✓ Branch 4 taken 170655 times.

      196460
                      if(d > dist[iv]) {
    
        1/2✓ Branch 1 taken 25805 times.
✗ Branch 2 not taken.

      25805
                          dist[iv] = d;
    
        1/2✓ Branch 1 taken 25805 times.
✗ Branch 2 not taken.

      25805
                          negative_pole[iv] = t;
    
                      }
    
                  }
    
              }
    
        4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.

      1
              Logger::out("LFS") << "(4) Kd-tree" << std::endl;
    
              // Step 4: create search structure
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.

      1
              vector<bool> is_pole(delaunay->nb_cells(), false);
    
        3/4✓ Branch 1 taken 10001 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10000 times.
✓ Branch 5 taken 1 times.

      10001
              for(index_t iv = 0; iv < delaunay->nb_vertices(); iv++) {
    
        3/4✓ Branch 1 taken 10000 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8476 times.
✓ Branch 4 taken 1524 times.

      10000
                  if(negative_pole[iv] != NO_POLE) {
    
        1/2✓ Branch 1 taken 8476 times.
✗ Branch 2 not taken.

      8476
                      is_pole[index_t(negative_pole[iv])] = true;
    
                  }
    
        2/2✓ Branch 2 taken 7105 times.
✓ Branch 3 taken 2895 times.

      10000
                  if(!voronoi_cell_is_infinite[iv]) {
    
        3/4✓ Branch 1 taken 7105 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6897 times.
✓ Branch 4 taken 208 times.

      7105
                      if(positive_pole[iv] != NO_POLE) {
    
        1/2✓ Branch 1 taken 6897 times.
✗ Branch 2 not taken.

      6897
                          is_pole[index_t(positive_pole[iv])] = true;
    
                      }
    
                  }
    
              }
    
      1
              index_t nb_poles = 0;
    
        3/4✓ Branch 1 taken 55452 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55451 times.
✓ Branch 5 taken 1 times.

      55452
              for(index_t t = 0; t < delaunay->nb_cells(); t++) {
    
        2/2✓ Branch 2 taken 15154 times.
✓ Branch 3 taken 40297 times.

      55451
                  if(is_pole[t]) {
    
      15154
                      nb_poles++;
    
                  }
    
              }
    
        1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.

      1
              poles_.reserve(nb_poles * 3);
    
        3/4✓ Branch 1 taken 55452 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 55451 times.
✓ Branch 5 taken 1 times.

      55452
              for(index_t t = 0; t < delaunay->nb_cells(); t++) {
    
        2/2✓ Branch 2 taken 15154 times.
✓ Branch 3 taken 40297 times.

      55451
                  if(is_pole[t]) {
    
        2/4✓ Branch 1 taken 15154 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15154 times.
✗ Branch 5 not taken.

      15154
                      poles_.push_back(circumcenter[t].x);
    
        2/4✓ Branch 1 taken 15154 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15154 times.
✗ Branch 5 not taken.

      15154
                      poles_.push_back(circumcenter[t].y);
    
        2/4✓ Branch 1 taken 15154 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 15154 times.
✗ Branch 5 not taken.

      15154
                      poles_.push_back(circumcenter[t].z);
    
                  }
    
              }
    
        3/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.

      1
              spatial_search_ = Delaunay::create(3, "NN");
    
        2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.

      1
              spatial_search_->set_vertices(poles_.size() / 3, poles_.data());
    
        4/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.

      1
              Logger::out("LFS") << "Done init." << std::endl;
    
      1
          }
    
      }

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#include <geogram/delaunay/LFS.h>
41			#include <geogram/basic/geometry.h>
42			#include <geogram/basic/logger.h>
43			#include <geogram/basic/command_line.h>
44
45			namespace {
46
47			using namespace GEO;
48
49			const index_t pid_size = 6;
50			int pid[pid_size][4] = {
51			{0, 1, 2, 3},
52			{0, 2, 1, 3},
53			{0, 3, 1, 2},
54			{1, 2, 0, 3},
55			{1, 3, 2, 0},
56			{2, 3, 0, 1}
57			};
58
59			/**
60			* \brief Computes the circumcenter and squared radius of a tetrahedron.
61			* \param[in] p first vertex of the tetrahedron
62			* \param[in] q second vertex of the tetrahedron
63			* \param[in] r third vertex of the tetrahedron
64			* \param[in] s fourth vertex of the tetrahedron
65			* \param[out] t_circumcenter computed circumcenter of the tetrahedron
66			* \param[out] squared_radius computed squared radius of the tetrahedron
67			* \param[out] dihedral_angles if non-nullptr,
68			* computed dihedral angles of the tetrahedron
69			* \return true if successful, false otherwise
70			* (for instance, if the tetrahedron is flat).
71			*/
72		55451	bool tetra_circumcenter_squaredradius(
73			const vec3& p, const vec3& q,
74			const vec3& r, const vec3& s,
75			vec3& t_circumcenter, double& squared_radius,
76			double* dihedral_angles = nullptr
77			) {
78
79		55451	vec3 qp = q - p;
80		55451	vec3 rp = r - p;
81		55451	vec3 sp = s - p;
82
83		55451	double tet_vol_x_12 = 2.0 * dot(qp, cross(rp, sp));
84
85			// it is not a safe check, one should avoid
86			// to input a degenerated tetra.
87	1/2 ✗ Branch 0 not taken. ✓ Branch 1 taken 55451 times.	55451	if(tet_vol_x_12 == 0.0) {
88		✗	return false;
89			}
90
91			vec3 d_t_circumcenter =
92	1/2 ✓ Branch 2 taken 55451 times. ✗ Branch 3 not taken.	55451	length2(qp) * cross(rp, sp) +
93	1/2 ✓ Branch 2 taken 55451 times. ✗ Branch 3 not taken.	110902	length2(rp) * cross(sp, qp) +
94	1/2 ✓ Branch 2 taken 55451 times. ✗ Branch 3 not taken.	110902	length2(sp) * cross(qp, rp);
95
96	1/2 ✓ Branch 1 taken 55451 times. ✗ Branch 2 not taken.	55451	squared_radius = length2(
97		55451	d_t_circumcenter) / (tet_vol_x_12 * tet_vol_x_12
98			);
99		55451	d_t_circumcenter = (1.0 / tet_vol_x_12) * d_t_circumcenter;
100		55451	t_circumcenter = p + d_t_circumcenter;
101
102	1/2 ✓ Branch 0 taken 55451 times. ✗ Branch 1 not taken.	55451	if(dihedral_angles) {
103			// compute Dihedral angle directly
104			// the following computation is not optimized.
105	2/2 ✓ Branch 1 taken 221804 times. ✓ Branch 2 taken 55451 times.	277255	vec3 point[4];
106		55451	point[0] = p;
107		55451	point[1] = q;
108		55451	point[2] = r;
109		55451	point[3] = s;
110	2/2 ✓ Branch 0 taken 332706 times. ✓ Branch 1 taken 55451 times.	388157	for(index_t i = 0; i < pid_size; i++) {
111		332706	vec3 b1 = point[pid[i][0]] - point[pid[i][2]];
112		332706	vec3 b2 = point[pid[i][1]] - point[pid[i][0]];
113		332706	vec3 b3 = point[pid[i][3]] - point[pid[i][1]];
114		332706	vec3 b2b3 = cross(b2, b3);
115		332706	dihedral_angles[i] = ::fabs(
116		332706	::atan2(
117	1/2 ✓ Branch 1 taken 332706 times. ✗ Branch 2 not taken.	332706	length(b2) * dot(b1, b2b3),
118		665412	dot(cross(b1, b2), b2b3)
119			)
120			);
121			}
122			}
123		55451	return true;
124			}
125
126			/**
127			* \brief Gets a Delaunay vertex by global index.
128			* \param[in] delaunay the Delaunay triangulation
129			* \param[in] v the index of the vertex
130			* \return a const reference to the vertex, as a vec3
131			*/
132		633474	inline const vec3& delaunay_vertex(Delaunay* delaunay, index_t v) {
133		633474	return (const vec3) delaunay->vertex_ptr(v);
134			}
135
136			/**
137			* \brief Gets a Delaunay vertex by tetrahedron index
138			* and local vertex index.
139			* \param[in] delaunay the Delaunay triangulation
140			* \param[in] c the index of the tetrahedron
141			* \param[in] lv the local index of the vertex (0,1,2 or 3)
142			* in tetrahedron \p c.
143			* \return a const reference to the vertex, as a vec3.
144			*/
145		221804	inline const vec3& delaunay_tet_vertex(
146			Delaunay* delaunay, index_t c, index_t lv
147			) {
148		221804	return delaunay_vertex(delaunay, index_t(delaunay->cell_vertex(c, lv)));
149			}
150
151			index_t tet_facet_vertex[4][3] = {
152			{1, 2, 3},
153			{0, 3, 2},
154			{3, 0, 1},
155			{2, 1, 0}
156			};
157
158			/**
159			* \brief Computes the normal to a tetrahedron facet.
160			* \param[in] delaunay the Delaunay triangulation
161			* \param[in] t the index of the tetrahedron
162			* \param[in] f the local index (0,1,2 or 3) of
163			* the facet in the tetrahedron \p t
164			* \return the normal to the facet \p f or tetrahedron \p t
165			*/
166		6250	vec3 delaunay_facet_normal(
167			Delaunay* delaunay, index_t t, index_t f
168			) {
169	1/6 ✗ Branch 0 not taken. ✓ Branch 1 taken 6250 times. ✗ Branch 3 not taken. ✗ Branch 4 not taken. ✗ Branch 6 not taken. ✗ Branch 7 not taken.	6250	geo_debug_assert(f < 4);
170		6250	index_t v1 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][0]));
171		6250	index_t v2 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][1]));
172		6250	index_t v3 = index_t(delaunay->cell_vertex(t, tet_facet_vertex[f][2]));
173		6250	const vec3& p1 = delaunay_vertex(delaunay, v1);
174		6250	const vec3& p2 = delaunay_vertex(delaunay, v2);
175		6250	const vec3& p3 = delaunay_vertex(delaunay, v3);
176		6250	return cross(p2 - p1, p3 - p1);
177			}
178			}
179
180			/****************************************************************************/
181
182			namespace GEO {
183
184		1	void LocalFeatureSize::init(index_t nb_pts, const double* pts) {
185			// Note: I need PDEL here instead of ANN/BNN since I need
186			// to access the cells (and ANN/BNN only give me the
187			// neighbors)
188
189	4/8 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 1 times. ✗ Branch 8 not taken. ✓ Branch 10 taken 1 times. ✗ Branch 11 not taken.	2	Logger::out("LFS") << "Delaunay" << std::endl;
190	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken.	1	Delaunay_var delaunay = Delaunay::create(3, "PDEL");
191	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken.	1	delaunay->set_vertices(nb_pts, pts);
192	4/8 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 1 times. ✗ Branch 8 not taken. ✓ Branch 10 taken 1 times. ✗ Branch 11 not taken.	1	Logger::out("LFS") << "Done Delaunay" << std::endl;
193
194	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 5 taken 1 times. ✗ Branch 6 not taken.	1	vector<vec3> circumcenter(delaunay->nb_cells());
195			vector<bool> voronoi_cell_is_infinite(
196			delaunay->nb_vertices(), false
197	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 5 taken 1 times. ✗ Branch 6 not taken.	1	);
198	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 5 taken 1 times. ✗ Branch 6 not taken.	1	vector<bool> is_sliver(delaunay->nb_cells());
199
200			static constexpr index_t NO_POLE = NO_INDEX;
201			// Index of incident tet whose circumcenter is the pole
202	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 5 taken 1 times. ✗ Branch 6 not taken.	1	vector<index_t> positive_pole(delaunay->nb_vertices(),NO_POLE);
203	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 5 taken 1 times. ✗ Branch 6 not taken.	1	vector<index_t> negative_pole(delaunay->nb_vertices(),NO_POLE);
204
205			vector<vec3> avg_infinite_dir(
206		1	delaunay->nb_vertices(), vec3(0, 0, 0)
207	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 5 taken 1 times. ✗ Branch 6 not taken.	1	);
208	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 5 taken 1 times. ✗ Branch 6 not taken.	1	vector<double> dist(delaunay->nb_vertices(), 0.0);
209
210	4/8 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 1 times. ✗ Branch 8 not taken. ✓ Branch 10 taken 1 times. ✗ Branch 11 not taken.	1	Logger::out("LFS") << "(1) Circumcenters and slivers" << std::endl;
211			// Step 1: compute circumcenters and check for slivers
212		1	const double sliver_quality = sliver_angle_threshold_ / 180.0 * M_PI;
213	3/4 ✓ Branch 1 taken 55452 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 55451 times. ✓ Branch 5 taken 1 times.	55452	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
214	1/2 ✓ Branch 2 taken 55451 times. ✗ Branch 3 not taken.	55451	const vec3& p = delaunay_tet_vertex(delaunay, t, 0);
215	1/2 ✓ Branch 2 taken 55451 times. ✗ Branch 3 not taken.	55451	const vec3& q = delaunay_tet_vertex(delaunay, t, 1);
216	1/2 ✓ Branch 2 taken 55451 times. ✗ Branch 3 not taken.	55451	const vec3& r = delaunay_tet_vertex(delaunay, t, 2);
217	1/2 ✓ Branch 2 taken 55451 times. ✗ Branch 3 not taken.	55451	const vec3& s = delaunay_tet_vertex(delaunay, t, 3);
218
219			double dihedral_angle[6];
220			double tet_radius;
221		110902	bool ok = tetra_circumcenter_squaredradius(
222			p, q, r, s,
223	2/4 ✓ Branch 1 taken 55451 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 55451 times. ✗ Branch 5 not taken.	55451	circumcenter[t], tet_radius, dihedral_angle
224			);
225
226	2/2 ✓ Branch 0 taken 332449 times. ✓ Branch 1 taken 55375 times.	387824	for(index_t a = 0; a < 6; a++) {
227	2/2 ✓ Branch 0 taken 76 times. ✓ Branch 1 taken 332373 times.	332449	if(!ok) {
228		76	break;
229			}
230		332373	ok = ok &&
231	3/4 ✓ Branch 0 taken 332373 times. ✗ Branch 1 not taken. ✓ Branch 2 taken 332295 times. ✓ Branch 3 taken 78 times.	664668	dihedral_angle[a] >= sliver_quality &&
232	2/2 ✓ Branch 0 taken 332287 times. ✓ Branch 1 taken 8 times.	332295	dihedral_angle[a] <= M_PI - sliver_quality;
233			}
234		55451	is_sliver[t] = !ok;
235			}
236
237	4/8 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 1 times. ✗ Branch 8 not taken. ✓ Branch 10 taken 1 times. ✗ Branch 11 not taken.	1	Logger::out("LFS") << "(2) Positive poles" << std::endl;
238			// Step 2: compute positive poles
239	3/4 ✓ Branch 1 taken 55452 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 55451 times. ✓ Branch 5 taken 1 times.	55452	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
240	2/2 ✓ Branch 2 taken 86 times. ✓ Branch 3 taken 55365 times.	55451	if(is_sliver[t]) {
241		86	continue;
242			}
243		55365	index_t f_inf = NO_INDEX;
244	2/2 ✓ Branch 0 taken 209269 times. ✓ Branch 1 taken 49115 times.	258384	for(index_t f = 0; f < 4; f++) {
245	4/6 ✓ Branch 1 taken 209269 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 209269 times. ✗ Branch 5 not taken. ✓ Branch 6 taken 6250 times. ✓ Branch 7 taken 203019 times.	209269	if(delaunay->cell_adjacent(t, f) == NO_INDEX) {
246		6250	f_inf = index_t(f);
247		6250	break;
248			}
249			}
250	2/2 ✓ Branch 0 taken 49115 times. ✓ Branch 1 taken 6250 times.	55365	if(f_inf == NO_INDEX) {
251			// tet t does not have facet on border
252	1/2 ✓ Branch 1 taken 49115 times. ✗ Branch 2 not taken.	49115	const vec3& c = circumcenter[t];
253	2/2 ✓ Branch 0 taken 196460 times. ✓ Branch 1 taken 49115 times.	245575	for(index_t lv = 0; lv < 4; lv++) {
254	2/4 ✓ Branch 1 taken 196460 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 196460 times. ✗ Branch 5 not taken.	196460	index_t iv = index_t(delaunay->cell_vertex(t, lv));
255	1/2 ✓ Branch 2 taken 196460 times. ✗ Branch 3 not taken.	196460	const vec3& v = delaunay_vertex(delaunay, iv);
256	1/2 ✓ Branch 2 taken 196460 times. ✗ Branch 3 not taken.	196460	double d = length2(c - v);
257	3/4 ✓ Branch 1 taken 196460 times. ✗ Branch 2 not taken. ✓ Branch 3 taken 34779 times. ✓ Branch 4 taken 161681 times.	196460	if(d > dist[iv]) {
258	1/2 ✓ Branch 1 taken 34779 times. ✗ Branch 2 not taken.	34779	dist[iv] = d;
259	1/2 ✓ Branch 1 taken 34779 times. ✗ Branch 2 not taken.	34779	positive_pole[iv] = t;
260			}
261			}
262			} else {
263			// tet t has facet f_inf on border
264			// positive pole = average direction of infinite Voronoi edges
265	2/4 ✓ Branch 1 taken 6250 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 6250 times. ✗ Branch 5 not taken.	6250	index_t v = index_t(delaunay->cell_vertex(t, index_t(f_inf)));
266		6250	voronoi_cell_is_infinite[v] = true;
267	1/2 ✓ Branch 1 taken 6250 times. ✗ Branch 2 not taken.	6250	avg_infinite_dir[v] +=
268	1/2 ✓ Branch 2 taken 6250 times. ✗ Branch 3 not taken.	12500	delaunay_facet_normal(delaunay, t, index_t(f_inf));
269			}
270			}
271
272	4/8 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 1 times. ✗ Branch 8 not taken. ✓ Branch 10 taken 1 times. ✗ Branch 11 not taken.	1	Logger::out("LFS") << "(3) Negative poles" << std::endl;
273			// Step 3: compute negative poles
274		1	std::fill(dist.begin(), dist.end(), 0.0);
275	3/4 ✓ Branch 1 taken 55452 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 55451 times. ✓ Branch 5 taken 1 times.	55452	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
276	2/2 ✓ Branch 2 taken 86 times. ✓ Branch 3 taken 55365 times.	55451	if(is_sliver[t]) {
277		86	continue;
278			}
279		55365	bool t_is_infinite = false;
280	2/2 ✓ Branch 0 taken 209269 times. ✓ Branch 1 taken 49115 times.	258384	for(index_t f = 0; f < 4; f++) {
281	4/6 ✓ Branch 1 taken 209269 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 209269 times. ✗ Branch 5 not taken. ✓ Branch 6 taken 6250 times. ✓ Branch 7 taken 203019 times.	209269	if(delaunay->cell_adjacent(t, f) == NO_INDEX) {
282		6250	t_is_infinite = true;
283		6250	break;
284			}
285			}
286	2/2 ✓ Branch 0 taken 6250 times. ✓ Branch 1 taken 49115 times.	55365	if(t_is_infinite) {
287		6250	continue;
288			}
289	1/2 ✓ Branch 1 taken 49115 times. ✗ Branch 2 not taken.	49115	const vec3& c = circumcenter[t];
290	2/2 ✓ Branch 0 taken 196460 times. ✓ Branch 1 taken 49115 times.	245575	for(index_t lv = 0; lv < 4; lv++) {
291	2/4 ✓ Branch 1 taken 196460 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 196460 times. ✗ Branch 5 not taken.	196460	index_t iv = index_t(delaunay->cell_vertex(t, lv));
292	1/2 ✓ Branch 2 taken 196460 times. ✗ Branch 3 not taken.	196460	const vec3& v = delaunay_vertex(delaunay, iv);
293		196460	vec3 N;
294	2/2 ✓ Branch 2 taken 40411 times. ✓ Branch 3 taken 156049 times.	196460	if(voronoi_cell_is_infinite[iv]) {
295	1/2 ✓ Branch 1 taken 40411 times. ✗ Branch 2 not taken.	40411	N = -avg_infinite_dir[iv];
296			} else {
297	2/4 ✓ Branch 1 taken 156049 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 156049 times. ✗ Branch 5 not taken.	156049	N = v - circumcenter[positive_pole[iv]];
298			}
299		196460	double d = dot(c - v, N);
300	3/4 ✓ Branch 1 taken 196460 times. ✗ Branch 2 not taken. ✓ Branch 3 taken 25805 times. ✓ Branch 4 taken 170655 times.	196460	if(d > dist[iv]) {
301	1/2 ✓ Branch 1 taken 25805 times. ✗ Branch 2 not taken.	25805	dist[iv] = d;
302	1/2 ✓ Branch 1 taken 25805 times. ✗ Branch 2 not taken.	25805	negative_pole[iv] = t;
303			}
304			}
305			}
306
307	4/8 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 1 times. ✗ Branch 8 not taken. ✓ Branch 10 taken 1 times. ✗ Branch 11 not taken.	1	Logger::out("LFS") << "(4) Kd-tree" << std::endl;
308			// Step 4: create search structure
309	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 5 taken 1 times. ✗ Branch 6 not taken.	1	vector<bool> is_pole(delaunay->nb_cells(), false);
310	3/4 ✓ Branch 1 taken 10001 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 10000 times. ✓ Branch 5 taken 1 times.	10001	for(index_t iv = 0; iv < delaunay->nb_vertices(); iv++) {
311	3/4 ✓ Branch 1 taken 10000 times. ✗ Branch 2 not taken. ✓ Branch 3 taken 8476 times. ✓ Branch 4 taken 1524 times.	10000	if(negative_pole[iv] != NO_POLE) {
312	1/2 ✓ Branch 1 taken 8476 times. ✗ Branch 2 not taken.	8476	is_pole[index_t(negative_pole[iv])] = true;
313			}
314	2/2 ✓ Branch 2 taken 7105 times. ✓ Branch 3 taken 2895 times.	10000	if(!voronoi_cell_is_infinite[iv]) {
315	3/4 ✓ Branch 1 taken 7105 times. ✗ Branch 2 not taken. ✓ Branch 3 taken 6897 times. ✓ Branch 4 taken 208 times.	7105	if(positive_pole[iv] != NO_POLE) {
316	1/2 ✓ Branch 1 taken 6897 times. ✗ Branch 2 not taken.	6897	is_pole[index_t(positive_pole[iv])] = true;
317			}
318			}
319			}
320
321		1	index_t nb_poles = 0;
322	3/4 ✓ Branch 1 taken 55452 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 55451 times. ✓ Branch 5 taken 1 times.	55452	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
323	2/2 ✓ Branch 2 taken 15154 times. ✓ Branch 3 taken 40297 times.	55451	if(is_pole[t]) {
324		15154	nb_poles++;
325			}
326			}
327
328	1/2 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken.	1	poles_.reserve(nb_poles * 3);
329	3/4 ✓ Branch 1 taken 55452 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 55451 times. ✓ Branch 5 taken 1 times.	55452	for(index_t t = 0; t < delaunay->nb_cells(); t++) {
330	2/2 ✓ Branch 2 taken 15154 times. ✓ Branch 3 taken 40297 times.	55451	if(is_pole[t]) {
331	2/4 ✓ Branch 1 taken 15154 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 15154 times. ✗ Branch 5 not taken.	15154	poles_.push_back(circumcenter[t].x);
332	2/4 ✓ Branch 1 taken 15154 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 15154 times. ✗ Branch 5 not taken.	15154	poles_.push_back(circumcenter[t].y);
333	2/4 ✓ Branch 1 taken 15154 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 15154 times. ✗ Branch 5 not taken.	15154	poles_.push_back(circumcenter[t].z);
334			}
335			}
336
337	3/6 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 1 times. ✗ Branch 8 not taken.	1	spatial_search_ = Delaunay::create(3, "NN");
338	2/4 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 6 taken 1 times. ✗ Branch 7 not taken.	1	spatial_search_->set_vertices(poles_.size() / 3, poles_.data());
339
340	4/8 ✓ Branch 1 taken 1 times. ✗ Branch 2 not taken. ✓ Branch 4 taken 1 times. ✗ Branch 5 not taken. ✓ Branch 7 taken 1 times. ✗ Branch 8 not taken. ✓ Branch 10 taken 1 times. ✗ Branch 11 not taken.	1	Logger::out("LFS") << "Done init." << std::endl;
341		1	}
342			}
343