minimal_shell_partitioner.h

/*
MIT License

Copyright (c) 2022 Michael Vennettilli

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
 */

#include <CGAL/Linear_cell_complex_for_combinatorial_map.h>
#include <CGAL/Polyhedron_3_to_lcc.h>
#include <CGAL/draw_linear_cell_complex.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Polyhedron_3.h>
#include <CGAL/convex_hull_3.h>

#include <vector>
#include <random>
#include <math.h>
#include <cmath>
#include <map>
#include <utility>

#ifndef MINIMAL_SHELL_PARTITIONER_H
#define MINIMAL_SHELL_PARTITIONER_H

// Typedefs pertaining to the kernel
typedef CGAL::Exact_predicates_inexact_constructions_kernel  K;
typedef CGAL::Polyhedron_3<K>                     Polyhedron_3;
typedef K::Point_3                                Point_3;
// Typedefs pertaining to the convex hull LCC
typedef CGAL::Linear_cell_complex_for_combinatorial_map<2,3> LCC_CH;
typedef LCC_CH::Dart_handle           Dart_handle_CH;
// Typedefs pertaining to the shell LCC
typedef CGAL::Linear_cell_complex_for_combinatorial_map<3> LCC_3;
typedef LCC_3::Dart_handle           Dart_handle_3;
typedef LCC_3::Vertex_attribute_handle           Vertex_handle_3;
// Typedefs pertaining to maps and pairs
typedef std::map<Dart_handle_CH, Vertex_handle_3>   CHDH_to_VH3_map;
typedef std::map<Dart_handle_CH, Dart_handle_3>   glue_vol_blueprint;
typedef std::pair<CHDH_to_VH3_map,CHDH_to_VH3_map>   CH_to_in_out_vertices;
typedef std::pair<std::vector<Dart_handle_3>, std::map<Dart_handle_3, Dart_handle_3>> inner_face_and_map;

inline int positive_modulo(int i, int base){
  return (base + (i % base)) % base;
}

/*
This function iterates over the vertex attributes for the convex hull
and prints out the points.
*/
void print_vertices(LCC_CH &lcc){
  for (LCC_CH::Vertex_attribute_range::iterator
         vertex_it=lcc.vertex_attributes().begin(),
         vertex_it_end=lcc.vertex_attributes().end();
       vertex_it!=vertex_it_end; ++vertex_it)
  {
    std::cout<<"point: "<<lcc.point_of_vertex_attribute(vertex_it)<<std::endl;
  }
}

/*
This computes the circumcenter of the face containing a given dart.
*/
LCC_3::Point circumcenter(LCC_CH &chull, Dart_handle_CH &dh){
  // Get the points
  LCC_CH::Point a = chull.point(dh);
  LCC_CH::Point b = chull.point(chull.beta(dh,1));
  LCC_CH::Point c = chull.point(chull.beta(dh,0));
  // Compute the differences
  LCC_CH::Vector ba = b-a;
  LCC_CH::Vector ca = c-a;
  // Compute the cross products
  LCC_CH::Vector X = CGAL::cross_product(ba,ca);
  LCC_CH::Vector Xba = CGAL::cross_product(X,ba);
  LCC_CH::Vector Xca = CGAL::cross_product(X,ca);
  // Compute the circumcenter
  LCC_3::Point cc = a + (Xba*ca.squared_length()-Xca*ba.squared_length())/(2*X.squared_length());
  return cc;
}

std::vector<Point_3> random_spherical_points(int num_pts){
  std::vector<Point_3> points;
  // Set up the random number generation
  std::mt19937 mt(time(NULL));
  double azimuthal_angle, polar_angle;
  const double pi = 3.14159265;
  std::uniform_real_distribution<double> azimuthal_dist(-1.0,1.0);
  std::uniform_real_distribution<double> polar_dist(0.0,2*pi);

  for(int i=0; i<num_pts; i++){
	  azimuthal_angle = acos(azimuthal_dist(mt));
	  polar_angle = polar_dist(mt);
	  points.push_back(Point_3(sin(azimuthal_angle)*cos(polar_angle),
			  sin(azimuthal_angle)*sin(polar_angle), cos(azimuthal_angle)));
  }
  return points;
}

/*
Computes the convex hull as a linear cell complex oriented clockwise
*/
LCC_CH make_chull(std::vector<Point_3> &points){
  // Compute the convex hull and convert to a linear cell complex
  Polyhedron_3 chull_polyhedron;
  CGAL::convex_hull_3(points.begin(), points.end(), chull_polyhedron);
  LCC_CH chull_lcc;
  Dart_handle_CH dh=CGAL::import_from_polyhedron_3<LCC_CH,Polyhedron_3>
    (chull_lcc, chull_polyhedron);
  // Check the orientation by computing the dot product between a cross product
  // generated by the ordered vertices of the face and the barycenter of that
  // face. Reverse the orientation if the dot product is positive.
  LCC_CH::Vector v1 = chull_lcc.point(chull_lcc.beta(dh,1))-chull_lcc.point(dh);
  LCC_CH::Vector v2 = chull_lcc.point(chull_lcc.beta(dh,0))-chull_lcc.point(dh);
  LCC_CH::Vector barycenter = chull_lcc.barycenter<2>(dh) - CGAL::ORIGIN;
  double signature = CGAL::cross_product(v1,v2)*barycenter;
  if(signature>0.0) chull_lcc.reverse_orientation();
  return chull_lcc;
}

/*
Use duality to get the vertices of the shell from the convex hull and a pair of
maps that send darts on a face of the convex hull to the corresponding inner/
outer vertex of the shell.
*/
CH_to_in_out_vertices make_dual_vertices(LCC_CH &chull,
  LCC_3 &shell, double r_in, double r_out){
  CHDH_to_VH3_map chull_to_inner_vertex;
  CHDH_to_VH3_map chull_to_outer_vertex;
  // Iterate over one dart per face in chull
  for(LCC_CH::One_dart_per_cell_range<2>::iterator
    face_it=chull.one_dart_per_cell<2>().begin(), face_it_end=chull.one_dart_per_cell<2>().end();
    face_it!=face_it_end;++face_it){
      // Compute the projections of the dual point onto the inner and outer spheres
      // Can use the barycenter of the face as an alternative
      LCC_3::Vector center = circumcenter(chull,face_it) - CGAL::ORIGIN;
      LCC_3::Point inner_point = CGAL::ORIGIN + (center * (r_in/std::sqrt(center.squared_length())));
      LCC_3::Point outer_point = CGAL::ORIGIN + (center * (r_out/std::sqrt(center.squared_length())));
      // Add the points to the vertex container and get their handles.
      Vertex_handle_3 inner_handle = shell.create_vertex_attribute<>(inner_point);
      Vertex_handle_3 outer_handle = shell.create_vertex_attribute<>(outer_point);
      // Iterate over the face and add the (dart, dual vertex) pairs to the
      // appropriate maps.
      Dart_handle_CH dh = face_it;
      do {
        chull_to_inner_vertex[dh] = inner_handle;
        chull_to_outer_vertex[dh] = outer_handle;
        dh = chull.beta(dh,1);
      } while (dh != face_it);
  }
  // Make the pair of maps and export it.
  CH_to_in_out_vertices map_pair = make_pair(chull_to_inner_vertex, chull_to_outer_vertex);
  return map_pair;
}

/*
This makes the inner and outer faces from one dart incident to a vertex in the
convex hull and updates a map that details how to glue the volumes. It returns
a vector containing the darts in the inner face and a map from darts in the
inner face to the outer face.
*/
inner_face_and_map make_inner_outer_pair(LCC_CH &chull, Dart_handle_CH &dh_ch_start,
  LCC_3 &shell, CH_to_in_out_vertices &map_pair, glue_vol_blueprint &glue_vol_map){
  std::vector<Dart_handle_3> inner_face, outer_face;
  std::map<Dart_handle_3, Dart_handle_3> inner_outer_map;
  Vertex_handle_3 inner_vh, outer_vh;
  Dart_handle_3 inner_dh, outer_dh;
  // Create a dart handle that will be used to iterate over the faces.
  Dart_handle_CH dh_ch = dh_ch_start;
  // Advance the faces until you've returned to the original one.
  do {
    // Get the vertex handles from the convex hull dart handle
    inner_vh = map_pair.first[dh_ch];
    outer_vh = map_pair.second[dh_ch];
    inner_dh = shell.create_dart(inner_vh);
    inner_face.push_back(inner_dh);
    outer_dh = shell.create_dart(outer_vh);
    outer_face.push_back(outer_dh);
    // Add the appropriate entry to glue_vol_map
    glue_vol_map[chull.beta(dh_ch,0)]=inner_dh;
    // Get the dart incident to the same vertex at the next (clockwise) face
    dh_ch = chull.beta(dh_ch, 0, 2);
  } while (dh_ch != dh_ch_start);
  // We make the association between the inner and outer darts and sew the faces.
  int vertices = inner_face.size();
  for(int i=0; i<vertices;++i){
    inner_outer_map[inner_face[i]] = outer_face[positive_modulo(i+1,vertices)];
    // The inner face should go clockwise
    shell.sew<1>(inner_face[i], inner_face[positive_modulo(i+1,vertices)]);
    // The outer face should go counterclockwise
    shell.sew<1>(outer_face[i], outer_face[positive_modulo(i-1,vertices)]);
  }
  inner_face_and_map result = make_pair(inner_face, inner_outer_map);
  return result;
}

/*
This function adds the lateral faces to connect the inner and outer faces
and close the volumes.
*/
void make_lateral_and_close(LCC_3 &shell, inner_face_and_map &inner_and_map_pair){
  std::vector<Dart_handle_3> inner_face = inner_and_map_pair.first;
  std::map<Dart_handle_3, Dart_handle_3> inner_outer_map = inner_and_map_pair.second;
  Dart_handle_3 l1, l2, l3, l4;
  Vertex_handle_3 vh1, vh2, vh3, vh4;
  int vertices = inner_face.size();
  for(int i=0; i<vertices; i++){
    // The dart 2-sewn to the inner dart starts at its destination.
    vh1 = shell.vertex_attribute(shell.beta(inner_face[i],1));
    l1 = shell.create_dart(vh1);
    // The dart after the new 2-sewn one has the same origin as the inner dart.
    vh2 = shell.vertex_attribute(inner_face[i]);
    l2 = shell.create_dart(vh2);
    // The dart 2-sewn to the outer dart starts at its destination.
    vh3 = shell.vertex_attribute(shell.beta(inner_outer_map[inner_face[i]],1));
    l3 = shell.create_dart(vh3);
    // The dart after the new 2-sewn one has the same origin as the outer dart.
    vh4 = shell.vertex_attribute(inner_outer_map[inner_face[i]]);
    l4 = shell.create_dart(vh4);
    // Connect the lateral darts to the inner and outer darts and make the face.
    shell.sew<2>(inner_face[i], l1);
    shell.sew<2>(inner_outer_map[inner_face[i]],l3);
    shell.sew<1>(l1,l2);
    shell.sew<1>(l2,l3);
    shell.sew<1>(l3,l4);
    shell.sew<1>(l4,l1);
  }
  // Now we 2-sew the lateral faces together.
  for(int i=0; i<vertices; i++){
    l1 = shell.beta(inner_face[i],2,1);
    l2 = shell.beta(inner_face[i],0,2,0);
    shell.sew<2>(l1,l2);
  }
  return;
}

/*
This is to be used after the closed volumes have been formed. This 3-sews
them together.
*/
void glue_vols(LCC_CH &chull, LCC_3 &shell, glue_vol_blueprint &glue_vol_map){
  Dart_handle_CH edge_it_flipped;
  Dart_handle_3 dh_sh1, dh_sh2;
  // Iterate over one dart per edge in the convex hull.
  for(LCC_CH::One_dart_per_cell_range<1>::iterator
    edge_it=chull.one_dart_per_cell<1>().begin(), edge_it_end=chull.one_dart_per_cell<1>().end();
    edge_it!=edge_it_end;++edge_it){
      // Get the other dart corresponding to the edge in the convex hull.
      edge_it_flipped = chull.beta(edge_it,2);
      // The glue_vol_map maps edges onto inner darts. To get the lateral
      // darts, you need to use beta_2.
      dh_sh1 = shell.beta(glue_vol_map[edge_it],2);
      dh_sh2 = shell.beta(glue_vol_map[edge_it_flipped],2);
      // 3-sew the darts together.
      shell.sew<3>(dh_sh1,dh_sh2);
  }
  return;
}

/*
This function puts everything together and uses duality to convert the convex
hull into a partition of the spherical shell.
*/
LCC_3 generate_shell(std::vector<Point_3> &points, double r_in, double r_out){
  LCC_3 shell;
  glue_vol_blueprint glue_vol_map;
  LCC_CH chull = make_chull(points);
  CH_to_in_out_vertices chd_to_shv = make_dual_vertices(chull, shell, r_in, r_out);
  // Iterate over one dart per vertex in chull
  for(LCC_CH::One_dart_per_cell_range<0>::iterator
    dart_it=chull.one_dart_per_cell<0>().begin(), dart_it_end=chull.one_dart_per_cell<0>().end();
    dart_it!=dart_it_end;++dart_it){
      // Make the inner and outer faces
      inner_face_and_map inner_and_map_pair = make_inner_outer_pair(chull, dart_it, shell,
        chd_to_shv, glue_vol_map);
      // Make the lateral edges and close the volume.
      make_lateral_and_close(shell, inner_and_map_pair);
    }
    // Glue the volumes and return the result.
    glue_vols(chull, shell, glue_vol_map);
    return shell;
}

/*
This computes the center of mass of a face.
*/
LCC_3::Point face_center_of_mass(LCC_3 &lcc, Dart_handle_3 dh_start){
  LCC_3::Point a, b;
  double weight = 0.0;
  LCC_3::Vector midpoint = LCC_3::Vector(0.0,0.0,0.0);
  LCC_3::Vector cumulative_sum = LCC_3::Vector(0.0,0.0,0.0);
  double perimeter = 0.0;
  Dart_handle_3 dh = dh_start;
  do {
    a = lcc.point(lcc.beta(dh,1));
    b = lcc.point(dh);
    weight = std::sqrt((b-a).squared_length());
    midpoint = ((a-CGAL::ORIGIN)+(b-CGAL::ORIGIN))*0.5;
    cumulative_sum += midpoint*weight;
    perimeter += weight;
    dh = lcc.beta(dh,1);
  } while(dh != dh_start);
  LCC_3::Point center_of_mass = CGAL::ORIGIN + (cumulative_sum/perimeter);
  return center_of_mass;
}

/*
This triangulates all faces in a LCC_3 by inserting the center of mass. Caution
needs to be used. Inserting a point into a face creates more faces, so use a mark
to indicate all new faces that have already been triangulated.
*/
void triangulate_all_faces(LCC_3 &lcc){
  LCC_3::size_type changed = lcc.get_new_mark();
  Dart_handle_3 dh_start, dh;
  LCC_3::Point center_of_mass;
  // Loop over all darts until you reach the (changing) end.
  LCC_3::Dart_range::iterator dart_it=lcc.darts().begin();
  while(dart_it!=lcc.darts().end()){
    // Have we already visited this dart?
    if(!lcc.is_marked(dart_it,changed)){
      // It belongs to a face we haven't divided yet. Split it.
      center_of_mass = face_center_of_mass(lcc,dart_it);
      dh_start = lcc.insert_point_in_cell<2>(dart_it, center_of_mass);
      dh = dh_start;
      // Mark all of the resulting faces belonging to the same volume as
      // dh_start
      while(!lcc.is_marked(dh,changed)){
        lcc.mark(dh,changed);
        lcc.mark(lcc.beta(dh,1),changed);
        lcc.mark(lcc.beta(dh,0),changed);
        dh = lcc.beta(dh,2,1);
      }
      // If dh_start is 3-sewn, we need to mark darts on the same face on
      // the other volume.
      if(!lcc.is_free<3>(dh_start)){
        dh = lcc.beta(dh_start,3,1);
        while(!lcc.is_marked(dh,changed)){
          lcc.mark(dh,changed);
          lcc.mark(lcc.beta(dh,1),changed);
          lcc.mark(lcc.beta(dh,0),changed);
          dh = lcc.beta(dh,2,1);
        }
      }
    }
    ++dart_it;
  }
  // Unmark everything before returning!
  lcc.unmark_all(changed);
  lcc.free_mark(changed);
  return;
}

/*
Compute a diagonal element for the matrix that arises in the generalization
of Lloyd's method to volumes on the sphere.
*/
double get_lloyd_scalar(LCC_3 &lcc, Dart_handle_3 dh, int i){
  // Get the ordered points, then compute the result.
  LCC_3::Point v0 = lcc.point(dh);
  LCC_3::Point v1 = lcc.point(lcc.beta(dh,1));
  LCC_3::Point v2 = lcc.point(lcc.beta(dh,0));
  double result = (v0[i]+v1[i]+v2[i])*(v0[i]+v1[i]+v2[i])
    - (v0[i]*v1[i]+v0[i]*v2[i]+v1[i]*v2[i]);
  return result;
}

/*
Given a dart on a face incident to a volume, this computes the diagonal matrix
and the outer normal, multiplies them together and returns the result.
*/
LCC_3::Vector lloyd_single_face(LCC_3 &lcc, Dart_handle_3 dh){
  LCC_3::Point v0 = lcc.point(dh);
  LCC_3::Point v1 = lcc.point(lcc.beta(dh,1));
  LCC_3::Point v2 = lcc.point(lcc.beta(dh,0));
  LCC_3::Vector outer_normal = CGAL::cross_product(v1-v0,v2-v0);
  double c0 = get_lloyd_scalar(lcc, dh, 0);
  double c1 = get_lloyd_scalar(lcc, dh, 1);
  double c2 = get_lloyd_scalar(lcc, dh, 2);
  LCC_3::Vector contribution = LCC_3::Vector(c0*outer_normal[0],
    c1*outer_normal[1],c2*outer_normal[2]);
  return contribution;
}

/*
Iterates over each volume in a 3-map and then over each face incident to that volume.
Computes the sum of the contributions from each face, normalizes the result, then
add export them as seeds for computing the convex hull.
*/
std::vector<LCC_CH::Point> get_seeds_lloyd(LCC_3 &lcc){
  std::vector<LCC_CH::Point> new_seeds;
  LCC_CH::Point seed;
  LCC_CH::Vector lloyd_sum;
  for(LCC_3::One_dart_per_cell_range<3>::iterator
    vol_it=lcc.one_dart_per_cell<3>().begin(), vol_end=lcc.one_dart_per_cell<3>().end();
    vol_it!=vol_end;++vol_it){
      // Reset the current seed and the sum.
      seed = LCC_CH::Point(0.0,0.0,0.0);
      lloyd_sum = LCC_CH::Vector(0.0,0.0,0.0);
      // Iterate over faces
      for(LCC_3::One_dart_per_incident_cell_range<2,3>::iterator
        face_it=lcc.one_dart_per_incident_cell<2,3>(vol_it).begin(),
        face_end=lcc.one_dart_per_incident_cell<2,3>(vol_it).end();
        face_it!=face_end; ++face_it){
          // Add the contribution from each face to the running sum
          lloyd_sum += lloyd_single_face(lcc,face_it);
        }
      // Normalize the sum, convert it to a point, then store it.
      lloyd_sum = lloyd_sum/std::sqrt(lloyd_sum.squared_length());
      seed = CGAL::ORIGIN + lloyd_sum;
      new_seeds.push_back(seed);
    }
    return new_seeds;
}

/*
This method recomputes the LCC_3 using Lloyd relaxation for a specified
number of iterations.
*/
void lloyd_relaxation(LCC_3 &lcc, int num_iter, double r_in, double r_out){
  std::vector<LCC_CH::Point> seeds;
  for(int i=0; i<num_iter; ++i){
    triangulate_all_faces(lcc);
    seeds = get_seeds_lloyd(lcc);
    lcc = generate_shell(seeds, r_in, r_out);
  }
  return;
}

#endif