Batch_collapser.h

Go to the documentation of this file.

// Copyright (c) 2012-2022 University of Lyon and CNRS (France).
// All rights reserved.
//
// This file is part of MEPP2; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 3 of
// the License, or (at your option) any later version.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 
#pragma once
 
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/properties.hpp>
 
#include "FEVV/Filters/Generic/generic_writer.hpp"
 
#include "CGAL/boost/graph/Euler_operations.h"
 
#include "FEVV/Filters/CGAL/Progressive_Compression/Metrics/Edge_length_metric.h"
#include "FEVV/Filters/CGAL/Progressive_Compression/Metrics/Volume_preserving.h"
#include "FEVV/Filters/CGAL/Progressive_Compression/Metrics/QEM_3D.h"
 
#include "FEVV/Tools/Comparator/Spanning_tree_vertex_edge_comparator.hpp"
 
#include "FEVV/Filters/CGAL/Progressive_Compression/Compression/Memory_comparator.h"
#include "FEVV/Filters/CGAL/Progressive_Compression/Predictors/Raw_positions.h"
#include "FEVV/Filters/CGAL/Progressive_Compression/Predictors/Butterfly.h"
 
#include "FEVV/Filters/CGAL/Progressive_Compression/apply_color.h"
 
#include "FEVV/Filters/CGAL/Progressive_Compression/Predictors/stencil.h"
#include "FEVV/Filters/CGAL/Progressive_Compression/Compression/Refinement_info.h"
#include "FEVV/Filters/CGAL/Progressive_Compression/Uniform_dequantization.h"
#include "FEVV/Filters/CGAL/Progressive_Compression/Helpers/Header_handler.h"
 
#include "FEVV/Filters/CGAL/Progressive_Compression/Geometric_metrics.h"
 
#include <iostream>
#include <fstream>
#include <vector>
#include <list>
#include <set>
#include <map>
#include <limits> 
 
namespace FEVV {
namespace Filters {
 
enum class BATCH_CONDITION { ALL_EDGES = 0, REACH_THRESHOLD };
 
template< typename HalfedgeGraph,
          typename PointMap,
          typename Metric,
          typename EdgeColorMap,
          typename VertexColorMap >
class Batch_collapser
{
public:
  using vertex_descriptor =
      typename boost::graph_traits< HalfedgeGraph >::vertex_descriptor;
  using halfedge_descriptor =
      typename boost::graph_traits< HalfedgeGraph >::halfedge_descriptor;
  using edge_descriptor =
      typename boost::graph_traits< HalfedgeGraph >::edge_descriptor;
  using face_descriptor =
      typename boost::graph_traits< HalfedgeGraph >::face_descriptor;
  using Vector = typename FEVV::Geometry_traits< HalfedgeGraph >::Vector;
  using Point = typename FEVV::Geometry_traits< HalfedgeGraph >::Point;
  using Geometry = typename FEVV::Geometry_traits< HalfedgeGraph >;
  
  Batch_collapser(
      HalfedgeGraph &g, 
      PointMap &pm, 
      Metric *metric, 
      Predictor< HalfedgeGraph, PointMap >
          *predictor,
      VertexColorMap &vcm, 
      EdgeColorMap &ecm, 
      BATCH_CONDITION batch_stop)
      : _g(g), _pm(pm), _metric(metric),
        _vcm(vcm), _ecm(ecm), 
        _gt(Geometry(g)), _predictor(predictor)
  {
    _batch_id = 0;
    _forbidden = 0;
    _link_condition = 0;
    _nb_collapse = 0;
    _number_vertices_original = FEVV::size_of_vertices(_g);
    _number_vertices_last = _number_vertices_original;
    _batch_stop = batch_stop;
  }
  ~Batch_collapser() {}
  
  bool stop()
  {
      switch(_batch_stop)
      {
      case BATCH_CONDITION::ALL_EDGES:
        return _metric->is_queue_empty();
        break;
      case BATCH_CONDITION::REACH_THRESHOLD:
        if(!_metric->is_queue_empty())
        {
 
          return (_metric->get_top_cost() > _metric->get_mean_threshold());
        }
        else
        {
          return true;
        }
        break;
 
      default:
        return _metric->is_queue_empty();
        break;
      }
 
  }
 
  void save_spanning_tree(
      const FEVV::Comparator::Spanning_tree_vertex_edge_comparator< HalfedgeGraph,
                                                          PointMap >
          &spanningtree)
  {
    std::ofstream file_encode;
    file_encode.open("encodingspanning" +
                     std::to_string(_batch_id) + ".txt");
    auto span = spanningtree.get_spanning_tree_vertices();
    typename std::list< vertex_descriptor >::const_iterator it = span.begin();
 
    for( ; it != span.end(); it++)
    {
      const Point& pt = get(_pm, *it);
      std::list< vertex_descriptor > list_v_around =
           FEVV::Comparator::get_adjacent_vertices(*it, _g);
      file_encode << _gt.get_x(pt) << " " << _gt.get_y(pt) << " "
                  << _gt.get_z(pt) << " " << list_v_around.size() << std::endl;
    }
    file_encode.close();
  }
 
  void compute_distortion_l2(
      FEVV::Filters::Geometric_metrics< HalfedgeGraph, PointMap > &g_metric, 
      Header_handler &header, 
      double diag, 
      bool skip 
                             )
  {
    if(!skip)
    {
      std::cout << "computing symmetrical L2 (max L2/Hausdorff and RMSE)" << std::endl;
      HalfedgeGraph current_graph = _g; // needs a copy 
      PointMap current_pm = get(boost::vertex_point, current_graph);
      Uniform_dequantization< HalfedgeGraph, PointMap > dq(
          current_graph,
          current_pm,
          header.get_quantization(),
          header.get_dimension(),
          header.get_init_coord());
 
      dq.point_dequantization();
      double geometric_error =
          g_metric.compute_symmetric_L2(current_graph, true);
      _RMSE_distortion_per_batch.push_back(geometric_error / diag);
      geometric_error =
          g_metric.compute_symmetric_L2(current_graph, false);
      _hausdorff_distortion_per_batch.push_back(geometric_error / diag);
    }
    else
    {
      _RMSE_distortion_per_batch.push_back(-1);
      _hausdorff_distortion_per_batch.push_back(-1);
    }
  }
 
  private:
  void sort_list_memory(
      const FEVV::Comparator::Spanning_tree_vertex_edge_comparator< HalfedgeGraph,
                                                          PointMap >
          &spanningtree)
  {
    Memory_comparator< HalfedgeGraph, PointMap > mc(spanningtree);
    _list_memory.sort(mc);
  }
  
    void setup_prediction(
      FEVV::Comparator::Spanning_tree_vertex_edge_comparator< HalfedgeGraph,
      PointMap >
    /*&spanningtree*/)
    {
      auto it_list = _list_memory.begin();
 
      for( ; it_list != _list_memory.end(); ++it_list)
      {
        _predictor->compute_residuals((*it_list));
      }
    }
 
public: 
  void collapse_batch()
  {
    // Count every case: forbidden edges, edges which would make our mesh non 
    // manifold or non-decodable.
    _forbidden = 0;
    _link_condition = 0;
    _metric->compute_error(); // Compute all edge costs for the current batch.
                              // Implements lines 7 to 11 of Algorithm 1.
 
    bool collapse_possible = true;
    _nb_collapse = 0;
 
    // Collapse until we meet our stopping condition.
    // Implements lines 12 to 23 of Algorithm 1.
    while(!stop())
    {
      collapse_possible = collapse_top_queue(); // Get the edge with lowest cost
                                                // and collapse it (if possible).
      if(collapse_possible)
      {
        _nb_collapse += 1;
      }
    }
    int64_t number_current_vertices = FEVV::size_of_vertices(_g);
#if(DEBUG)
    // Compute percentage of removed vertices compared to original mesh.
    double percentage_vertices = (1 - (double)number_current_vertices /
                                         (double)_number_vertices_original) * 100;
    std::cout << "Percentage vertices removed from start " << percentage_vertices << std::endl;
 
    percentage_vertices = (1 - (double)number_current_vertices / (double)_number_vertices_last) * 100;
    std::cout << "Percentage vertices removed from last " << percentage_vertices << std::endl;
#endif
    _number_vertices_last = number_current_vertices;
    // Spanning tree construction. Implements line 24 of Algorithm 1.
    FEVV::Comparator::Spanning_tree_vertex_edge_comparator< HalfedgeGraph,
                                                        PointMap > spanningtree(_g, _pm, true);
    sort_list_memory(spanningtree); // Sort _list_memory according to vertex 
                                    // st traversal.
                                    // Implements line 25 of Algorithm 1.
    // Store refinement info as bitstreams and residuals.
    Refinement_info< HalfedgeGraph,
                    PointMap >
        ref_settings(_g, _list_memory);
 
    // set_bitMask, set_connectivity_topology, and set_reverse_bool
    // implement line 26 of Algorithm 1.
    // Compute vertex to split bitmask. 
    ref_settings.set_bitMask(spanningtree); // set vertex to split bitmask
 
    // Compute connectivity bitmask (to know which halfedges to expand at the
    // decompression step).
    ref_settings.set_connectivity_topology(spanningtree, 
                                           _pm, 
                                           _list_memory); // set one-ring bitmask
    ref_settings.set_reverse_bool();// set reverse bitmask
 
    // Compute Predictions.
    // The two lines below implement line 27 of Algorithm 1.
    setup_prediction(spanningtree); // compute residuals
    ref_settings.set_error_prediction(); // set residuals array
    _refinements.push_back(ref_settings); // Implements line 28 of Algorithm 1
#if(DEBUG)
    std::cout << "batch nb " << _batch_id << " done" << std::endl;
#endif
    _batch_id++;
    _forbidden_edges.clear();
    _list_memory.clear();
  }
 
  private:
  bool are_halfedge_vertex_positions_change_without_normal_flip(halfedge_descriptor h_collapse,
                           const Point &pos_vkept)
  {
    return (is_vertex_position_change_without_normal_flip(h_collapse, source(h_collapse, _g), pos_vkept) &&
            is_vertex_position_change_without_normal_flip(h_collapse, target(h_collapse, _g), pos_vkept));
  }
 
  bool is_vertex_position_change_without_normal_flip(
                                                halfedge_descriptor h_collapse,
                                                vertex_descriptor vertex,
                                                const Point& pos_vkept)
  {
    assert(vertex == source(h_collapse, _g) || 
         vertex == target(h_collapse, _g));
 
    const Point& pvertex = get(_pm, vertex);
 
    if(pos_vkept == pvertex)
      return true;
    else
    {
      auto iterator_range = CGAL::halfedges_around_target(vertex, _g);
      for(auto h : iterator_range)
      {
        // case where h is a border
        if(CGAL::is_border(h, _g))
          continue;
 
        face_descriptor f = face(h, _g); // we know that f!=null_face()
 
        // when f is one of the two incident faces to the collapsed edge
        if(f == face(h_collapse, _g))
          continue;
 
        if(f == face(opposite(h_collapse, _g), _g))
          continue;
 
        const Point& p1 = get(_pm, source(h, _g));
        const Point& p0 = get(_pm, source(prev(h, _g), _g));
 
        if((p0 == pvertex) || (p1 == pvertex))
        {
          std::cerr << "Warning: is_vertex_position_change_without_normal_flip: "
                       "we should not have dupplicate at this stage ";
        }
 
        Vector face_normal, face_normal_after_collapse;
        // get the normal of the current face
        // check if the current face is flat
        bool b = FEVV::Math::Vector::are_collinear< Geometry >(
            _gt.sub_p(p0, p1), _gt.sub_p(pvertex, p1));
        if(!b)
          face_normal =
              _gt.normal(p0, p1, pvertex); // gt.unit_normal(p0, p1, pvertex);
        else                               
        { // when the current face is flat
          face_normal =
              Vector(std::numeric_limits< typename FEVV::Geometry_traits<
                         HalfedgeGraph >::Scalar >::quiet_NaN(),
                     std::numeric_limits< typename FEVV::Geometry_traits<
                         HalfedgeGraph >::Scalar >::quiet_NaN(),
                     std::numeric_limits< typename FEVV::Geometry_traits<
                         HalfedgeGraph >::Scalar >::quiet_NaN());
 
          //return false; // the face was flat, but with the new pos_vkept
                      // can resolve it
        }
        // simulation of the new face after collapse and check if its
        // normal has the same sign 
        auto cross = _gt.cross_product(p0 - p1, p0 - pos_vkept);
        if(cross != Vector(0, 0, 0))
        {
          face_normal_after_collapse = _gt.unit_normal(
              p0, p1, pos_vkept);
 
          const typename Geometry::Scalar sign =
              _gt.dot_product(face_normal, face_normal_after_collapse);
          if(sign < 0.0)
            return false;
        }
        else
        {
          return false;
        }
      }
 
      return true;
    }
  }
 
  bool is_forbidden(halfedge_descriptor h, Point /*pvkept*/)
  {
    // forbid border cases for the first version
    return (
        ((_forbidden_edges.find(h) == _forbidden_edges.end() &&
          _forbidden_edges.find(opposite(h, _g)) == _forbidden_edges.end()) &&
         _metric->is_present_in_map(edge(h, _g))));
  }
 
  bool collapse_top_queue()
  {
    bool mesh_changed = false;
    if(!_metric->is_queue_empty())
    {
      // we obtain the edge with the lowest cost
      std::tuple< edge_descriptor, double, Point > current_edge =
          _metric->get_top_queue(); // Implements line 13 of Algorithm 1.
      
      // is the  current edge forbidden?
      auto current_halfedge = halfedge(std::get< 0 >(current_edge), _g);
 
      if(is_forbidden(current_halfedge, std::get< 2 >(current_edge)))
      {
        if(CGAL::Euler::does_satisfy_link_condition(std::get< 0 >(current_edge),
                                                    _g))
        {
          vertex_descriptor vs = source(current_halfedge, _g);
          vertex_descriptor vt = target(current_halfedge, _g);
 
          const Point& pvkept = std::get< 2 >(current_edge); // vkept position
 
          // Check if collapsing this edge would flip normals, thus adding
          // artifacts. If it does, do not collapse.
          bool normal_flip = are_halfedge_vertex_positions_change_without_normal_flip(current_halfedge, pvkept); // true when ok
 
          if(normal_flip)
          {
            // Here the edge is collapsible (line 14 of Algorithm 1).
 
            // Create a collapse info objet to store the info we need to
            // reconstruct the edge.
            Collapse_info< HalfedgeGraph, PointMap > memory(_g, _pm);
            memory.record_vt_vs_pos(get(_pm, vt), get(_pm, vs)); 
            memory.record_v3_v4(current_halfedge); 
 
            // give collapse info a unique ID (useful for debugging)
            memory.set_num_collapse(_nb_collapse);
 
            // forbid the simplification of the neighbourhood of the collapsed
            // edge for the current batch
            forbid_edges(_g, vs, _forbidden_edges); // Implements line 16 of Algorithm 1.
            forbid_edges(_g, vt, _forbidden_edges); // Implements line 17 of Algorithm 1.           
            
            // Collapse edge.
            // Implements line 18 of Algorithm 1.
            vertex_descriptor vkept =
                  CGAL::Euler::collapse_edge(std::get< 0 >(current_edge), _g);
            mesh_changed = true;
 
            // set position of the kept vertex
            put(_pm, vkept, pvkept);
 
            memory.record_pos_vkept(vkept);
            memory.record_vkept(vkept);
 
            _list_memory.push_back(std::move(memory));
          }
#if(DEBUG)      
          else
          {
              std::cerr << "Warning: collapse_top_queue: an collapse did not occur due to a normal flip." << std::endl;
          }
#endif
          _link_condition++;
        }
      }
      else
      {
        _forbidden++;
        mesh_changed = false;
      }
      _metric->pop_queue(); // dans tous les cas on pop
    }
    return mesh_changed;
  }
 
public:
  const std::vector< Refinement_info< HalfedgeGraph,
                               PointMap > >&
  get_refinements() const
  {
 
    return _refinements;
  }
 
  const std::vector< double >& get_RMSE_distortion() const { return _RMSE_distortion_per_batch; }
  const std::vector< double >& get_hausdorff_distortion() const { return _hausdorff_distortion_per_batch; }
private:
  HalfedgeGraph &_g; 
  PointMap &_pm; 
  Metric *_metric; 
  VertexColorMap &_vcm;
  EdgeColorMap &_ecm;
  std::set< halfedge_descriptor > _forbidden_edges; 
  const Geometry _gt;
  int _batch_id; 
  int _nb_collapse;
  std::list< Collapse_info< HalfedgeGraph, PointMap > > _list_memory; 
  FEVV::Filters::
      Predictor< HalfedgeGraph, PointMap >
          *_predictor;
  std::vector< Refinement_info< HalfedgeGraph,
                               PointMap > >
      _refinements; 
 
  std::vector< double > _RMSE_distortion_per_batch, _hausdorff_distortion_per_batch;
  std::vector< std::vector< bool > > _list_edge_bitmask;
  int _link_condition;
  int _forbidden;
  int64_t _number_vertices_original, _number_vertices_last;
  BATCH_CONDITION _batch_stop;
};
 
} // namespace Filters
} // namespace FEVV