partition_optimal_convex_2.h

很多二维 三维几何计算算法 C++ 类库

// Copyright (c) 2000  Max-Planck-Institute Saarbruecken (Germany).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org); you may redistribute it under
// the terms of the Q Public License version 1.0.
// See the file LICENSE.QPL distributed with CGAL.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// $URL: svn+ssh://scm.gforge.inria.fr/svn/cgal/branches/CGAL-3.3-branch/Partition_2/include/CGAL/Partition_2/partition_optimal_convex_2.h $
// $Id: partition_optimal_convex_2.h 37992 2007-04-07 09:39:45Z efif $
// 
//
// Author(s)     : Susan Hert <[e-mail]>


// ===========================================================================
// There is a known bug in this algorithm that causes more than the optimal 
// number of convex pieces to be reported in cases where there are many
// collinear vertices (as in a Hilbert polygon, for example).  More precisely,
// the problem is known to crop up in this situation:
//
//       5-----------4
//       |           |
//       |     2-----3
//       |     |
//       |     1-----0
//       |           |
//       |    13-----14
//       |     |
//       6    12--11
//       |        |
//       |     9--10
//       |     |
//       7-----8
//
// The problem arises because when decomposing the polygon from vertex 1 to 13
// point 2 is (quite correctly) indicated as not visible to point 13.  Thus
// it is believed an edge is necessary to divide the polygon (1 2 5 6 12 13).
//
// A hack that partially fixes this problem is implemented as follows:
// a vertex r is marked as visible from a point q for the purposes of 
// the decompose function if p is the other endpoint of the edge containing r 
// and p is visible from q.
//
// This causes the problem that decomposition from 8 to 12 indicates that 
// valid vertices are 8 9 and 12. Diagonal (9 12) is valid, but (8 12) is
// also considered to be valid and necessary since 9 is a reflex vertex. 
// The result is a polygon split with diagonals (2 5) (12 5) (9 12) and
// (13 1), which is obviously not optimal.
//
// To get around this problem, I currently postprocess during the cutting
// up of the polygon to remove unneeded diagonals and then
// achieve the optimal, but a better solution is surely desired....
//
// ===========================================================================

#ifndef CGAL_PARTITION_OPTIMAL_CONVEX_H
#define CGAL_PARTITION_OPTIMAL_CONVEX_H

#include<CGAL/IO/Tee_for_output_iterator.h>
#include<CGAL/Partition_2/Partition_opt_cvx_edge.h>
#include<CGAL/Partition_2/Partition_opt_cvx_vertex.h>
#include<CGAL/Partition_2/Partition_opt_cvx_diagonal_list.h>
#include<CGAL/Partition_2/Matrix.h>
#include<CGAL/Partition_2/Turn_reverser.h>
#include<CGAL/Partition_2/Partitioned_polygon_2.h>
#include<CGAL/partition_is_valid_2.h>
#include<CGAL/Partition_traits_2.h>
#include<CGAL/Partition_2/partition_assertions.h>
#include<CGAL/Partition_2/Vertex_visibility_graph_2.h>
#include<utility>
#include<vector>
#include<iterator>
#include<iostream>


namespace CGAL {


#ifdef  CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
int partition_opt_cvx_debug_list_count = 0;  
#endif


template <class Polygon, class Traits>
int partition_opt_cvx_best_so_far(Partition_opt_cvx_vertex& pivot_vertex,
                                  unsigned int extension,
                                  Polygon& polygon, const Traits& traits,
                                  Partition_opt_cvx_diagonal_list& diag_list)
{
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   std::cout << "best(" << pivot_vertex.vertex_num() << "." 
             << partition_opt_cvx_debug_list_count 
             << ", " << extension << ")" << std::endl;
#endif
   Partition_opt_cvx_stack_record best_so_far = pivot_vertex.best_so_far();
   while (!pivot_vertex.stack_empty()) {
      Partition_opt_cvx_stack_record old = pivot_vertex.stack_top();
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
      std::cout << "best(" << pivot_vertex.vertex_num() << "." 
                << partition_opt_cvx_debug_list_count << ", " << extension 
                << ")" << " old = " << old.vertex_num() 
                << ", " << old.value() << std::endl;
#endif
      typedef typename Traits::Left_turn_2    Left_turn_2;
      typedef typename Traits::Point_2       Point_2;
      Left_turn_2 left_turn = traits.left_turn_2_object();
      Turn_reverser<Point_2, Left_turn_2>  right_turn(left_turn);
      if (right_turn(polygon[old.vertex_num()], 
                    polygon[pivot_vertex.vertex_num()], polygon[extension])) 
      {
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
         std::cout << "best(" << pivot_vertex.vertex_num() << "." 
                   << partition_opt_cvx_debug_list_count 
                   << ", " << extension << ")" << " returning(a) " 
                   << best_so_far.value() << std::endl;
#endif
         diag_list = best_so_far.solution();
         return best_so_far.value();
      }
      else if (old.value() < best_so_far.value())
         best_so_far = old;
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
      std::cout << "best(" << pivot_vertex.vertex_num() << "." 
                << partition_opt_cvx_debug_list_count 
                << ", " << extension << ") " << "popping off " 
                << old.vertex_num() << ", " << old.value() << std::endl;
#endif
      pivot_vertex.stack_pop();
   }
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   std::cout << "best(" << pivot_vertex.vertex_num() << "." 
             << partition_opt_cvx_debug_list_count 
             << ", " << extension << ") returning(b) " << best_so_far.value() 
             << std::endl;
#endif
   diag_list = best_so_far.solution();
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   std::cout << "     diagonal list " << diag_list << std::endl; 
#endif
   return best_so_far.value();
}

template <class Polygon, class Traits>
void partition_opt_cvx_load(int current, 
                            ::std::vector< Partition_opt_cvx_vertex >& v_list,
                            Polygon& polygon, 
                            Matrix<Partition_opt_cvx_edge>& edges, 
                            const Traits& traits)
{
    int previous;
    int num_polygons;
    Partition_opt_cvx_diagonal_list diag_list1, diag_list2;

#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
    std::cout << "load(" << v_list[current].vertex_num() << ")" << std::endl;
#endif
    for (previous = current-1; previous >= 0; previous--) {
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
      std::cout << "load:  previous = " << v_list[previous].vertex_num() 
                << std::endl;
#endif
      // must look at all valid edges and at all edges that are visible and
      // have something on the stack.  The latter check is necessary to make
      // sure solutions accumulate properly.
      if (edges[v_list[previous].vertex_num()]
               [v_list[current].vertex_num()].is_valid() ||
          (edges[v_list[previous].vertex_num()]
                [v_list[current].vertex_num()].is_visible() &&
           !v_list[previous].stack_empty())) 
      {
         num_polygons = partition_opt_cvx_decompose(
                                  v_list[previous].vertex_num(), 
                                  v_list[current].vertex_num(), polygon, 
                                  edges, traits, diag_list1) + 
                        partition_opt_cvx_best_so_far(v_list[previous], 
                                                  v_list[current].vertex_num(),
                                                  polygon, traits, diag_list2);
         diag_list1.splice(diag_list1.end(), diag_list2);
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
         std::cout << "load:  pushing previous = " 
                   << v_list[previous].vertex_num() << " num_polygons = " 
                   << num_polygons << " on stack " 
                   << v_list[current].vertex_num() << "." 
                   << partition_opt_cvx_debug_list_count << std::endl;
         std::cout << "     diagonal list = " << diag_list1 << std::endl;
#endif
         v_list[current].stack_push(v_list[previous].vertex_num(), 
                                    num_polygons, diag_list1);
       }
    }
}

// pre: edge_num1 <= e_num <= edge_num2 but edge_num1 != edge_num2
template <class Polygon, class Traits>
bool collinearly_visible(unsigned int edge_num1, unsigned int e_num,
                         unsigned int edge_num2,
                         const Matrix<Partition_opt_cvx_edge>& edges,
                         const Polygon& polygon,
                         const Traits& traits)
{
   typedef typename Polygon::size_type                   size_type;
   typedef typename Traits::Orientation_2                Orientation_2;
   Orientation_2 orientation = traits.orientation_2_object();

   if ((e_num == edge_num1+1 || e_num+1 == edge_num2) && 
       edges[edge_num1][edge_num2].is_visible() &&
       orientation(polygon[edge_num1], polygon[e_num], polygon[edge_num2]) ==
                  COLLINEAR)
     return true;
   else
     return false;    
}
           

template <class Polygon, class Traits>
int partition_opt_cvx_decompose(unsigned int edge_num1, unsigned int edge_num2,
                                Polygon& polygon, 
                                Matrix<Partition_opt_cvx_edge>& edges, 
                                const Traits& traits, 
                                Partition_opt_cvx_diagonal_list& diag_list)
{
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   std::cout << "decompose(" << edge_num1 << ", " << edge_num2 << ")";
#endif
   if (edges[edge_num1][edge_num2].is_done())  {
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
      std::cout << " returning " << edges[edge_num1][edge_num2].value() 
                << std::endl;
#endif
      diag_list = edges[edge_num1][edge_num2].solution();
      return edges[edge_num1][edge_num2].value();
   }
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   std::cout << std::endl;
#endif

   // temporarily invalidate this edge so we don't try to decompose on this
   // edge again 
   Partition_opt_cvx_edge_validity old_validity; 
   old_validity = edges[edge_num1][edge_num2].validity();
   edges[edge_num1][edge_num2].set_valid(PARTITION_OPT_CVX_NOT_VALID);

   std::vector< Partition_opt_cvx_vertex > v_list; 
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   partition_opt_cvx_debug_list_count++;
#endif
   typedef typename Polygon::size_type  size_type;

   for (size_type e_num = edge_num1; e_num <= edge_num2; e_num++) 
   {
       if ((edges[edge_num1][e_num].is_visible() && 
            edges[e_num][edge_num2].is_visible() ) || 
           collinearly_visible(edge_num1, e_num, edge_num2, edges, polygon,
                               traits) ) 
       {
          v_list.push_back(Partition_opt_cvx_vertex(e_num));
       }
   }
   std::vector< int >::size_type v;

#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   std::cout << "v_list(" << edge_num1 << ", " << edge_num2 << ")";
   for(v = 0; v < v_list.size(); v++) {
       std::cout << " " << v_list[v].vertex_num();
   }
   std::cout << std::endl;
#endif

   for(v = 0; v < v_list.size(); v++) {
       partition_opt_cvx_load(int(v), v_list, polygon, edges, traits);
   }

   int num_pieces = partition_opt_cvx_best_so_far(v_list[v_list.size()-1], 
                                                  edge_num1, polygon, traits,
                                                  diag_list) + 1;
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   std::cout << "decompose: num_pieces = " << num_pieces << std::endl;
#endif
   edges[edge_num1][edge_num2].set_value(num_pieces);
   diag_list.push_back(Partition_opt_cvx_diagonal(edge_num1, edge_num2));
   edges[edge_num1][edge_num2].set_value(num_pieces);
   edges[edge_num1][edge_num2].set_solution(diag_list);
   edges[edge_num1][edge_num2].set_done(true);
   edges[edge_num1][edge_num2].set_valid(old_validity);
   // revalidate the edge; next time it will pick up the computed value 
   // stored with this edge
#ifdef CGAL_PARTITION_OPTIMAL_CONVEX_DEBUG
   std::cout << "decompose(" << edge_num1 << ", " << edge_num2 << "): "
             << " edge[" << edge_num1 << "][" << edge_num2 << "] set to " 
             << edges[edge_num1][edge_num2] << std::endl;
   std::cout << " with diagonal list " 
             << edges[edge_num1][edge_num2].solution()