nesting.cpp

#include "nesting.h"

namespace nesting {
NFPCacheValue& comp_nfp(const Polygon_with_holes_2* poly_A,
                        const uint32_t rotation_A,
                        const uint32_t allowed_rotation_A,
                        const Polygon_with_holes_2* poly_B,
                        const uint32_t rotation_B,
                        const uint32_t allowed_rotation_B,
                        Layout& layout) {
  NFPCacheKey key(poly_A, poly_B, rotation_A, rotation_B);
  Polygon_with_holes_2 nfp;
  auto kv = layout.nfp_cache.find(key);
  if (kv == layout.nfp_cache.end()) {
    Transformation scale(CGAL::SCALING, -1);
    auto rotate_A = geo::get_rotate(rotation_A, allowed_rotation_A);
    auto rotate_B = geo::get_rotate(rotation_B, allowed_rotation_B);
    Polygon_with_holes_2 minus_B;
    Polygon_with_holes_2 rotated_A;
    if (rotate_B) {
      minus_B = geo::transform_polygon_with_holes(scale * (*rotate_B), *poly_B);
    } else {
      minus_B = geo::transform_polygon_with_holes(scale, *poly_B);
    }
    if (rotate_A) {
      rotated_A = geo::transform_polygon_with_holes(*rotate_A, *poly_A);
    } else {
      rotated_A = *poly_A;
    }
    nfp = CGAL::minkowski_sum_2(rotated_A, minus_B);
    geo::strict_simplify(nfp);
    auto bbox = nfp.bbox();
    NFPCacheValue v;
    v.nfp = nfp;
    v.xmin = CGAL::to_double(bbox.xmin());
    v.xmax = CGAL::to_double(bbox.xmax());
    v.ymin = CGAL::to_double(bbox.ymin());
    v.ymax = CGAL::to_double(bbox.ymax());
    auto _kv = layout.nfp_cache.emplace(key, v);
    return _kv.first->second;
  } else {
    return kv->second;
  }
}
FT comp_pd(NFPCacheValue& v, double px, double py, Layout& layout) {
  auto& original_nfp = v.nfp;
  if (px <= v.xmin || px >= v.xmax || py <= v.ymin || py >= v.ymax) {
    return geo::FT_ZERO;
  }
  hash::PDCacheKey key(&original_nfp, px, py);
  FT pd;
  layout.pd_count++;
  auto iter = layout.pd_cache.find(key);
  if (iter != layout.pd_cache.end()) {
    pd = iter->second;
  } else {
    layout.pd_miss++;
    // 为效率考虑使用不精确构建
    Point_2 relative_point(px, py);
    pd = geo::comp_pd(original_nfp, relative_point);
    layout.pd_cache.emplace(key, pd);
  }
  return pd;
}
void shrink(Layout& layout) {
  layout.cur_length = layout.best_length * (1 - layout.rdec);
  if (layout.cur_length < layout.lower_length) {
    layout.cur_length = layout.lower_length;
  }
  std::vector<size_t> v;
  // 将超出cur_length的polygon重新放置
  for (size_t i = 0; i < layout.poly_num; i++) {
    auto& p = layout.sheet_parts[0][i];
    auto right_vertex = p.transformed.outer_boundary().right_vertex();
    auto x = right_vertex->x();
    if (x > layout.cur_length) {
      auto ifr = geo::comp_ifr(layout.sheets[0].sheet, p.transformed);
      auto bbox = ifr.bbox();
      double random_x;
      double random_y;
      // 放置策略：随机放置
      random_x = bbox.x_span() * rand::right_nd01() + bbox.xmin();
      random_y = bbox.y_span() * rand::center_nd01() + bbox.ymin();
      p.set_translate(random_x, random_y);
      v.push_back(i);
    }
  }
  // 更新pd
  for (auto& i : v) {
    auto& p = layout.sheet_parts[0][i];
    auto& px = p.get_translate_ft_x();
    auto& py = p.get_translate_ft_y();
    auto double_px = CGAL::to_double(px);
    auto double_py = CGAL::to_double(py);
    for (size_t j = 0; j < layout.poly_num; j++) {
      if (i == j) {
        continue;
      }
      auto& q = layout.sheet_parts[0][j];
      auto& qx = q.get_translate_ft_x();
      auto& qy = q.get_translate_ft_y();
      auto double_qx = q.get_translate_double_x();
      auto double_qy = q.get_translate_double_y();
      auto& nfp =
          comp_nfp(q.base, q.get_rotation(), q.allowed_rotations, p.base,
                   p.get_rotation(), p.allowed_rotations, layout);
      Point_2 relative_point(px - qx, py - qy);
      auto pd =
          comp_pd(nfp, double_px - double_qx, double_py - double_qy, layout);
      layout.set_pd(i, j, pd);
    }
  }
  // 更新sheet
  layout.sheets[0].set_width(layout.cur_length);
}
void get_init_solu(Layout& layout) {
  //std::ofstream s("init.txt", std::ios::out);
  //s << "init solu start" << std::endl;
  size_t num_poly = layout.poly_num;
  for (size_t i = 0; i < num_poly; i++) {
    auto& shape_i = layout.sheet_parts[0][i];
    auto& polygon_i = shape_i.transformed;
    auto rotation_i = shape_i.get_rotation();
    auto& ix = shape_i.get_translate_ft_x();
    auto& iy = shape_i.get_translate_ft_y();
    auto double_ix = CGAL::to_double(ix);
    auto double_iy = CGAL::to_double(iy);
    auto ifr = geo::comp_ifr(layout.sheets[0].sheet, polygon_i);
    CandidatePoints c(layout.poly_num - 1);
    c.set_boundary(ifr);
    for (size_t j = 0; j < num_poly; j++) {
      if (j == i) {
        continue;
      }
      auto& shape_j = layout.sheet_parts[0][j];
      auto& jx = shape_j.get_translate_ft_x();
      auto& jy = shape_j.get_translate_ft_y();
      auto double_jx = shape_j.get_translate_double_x();
      auto double_jy = shape_j.get_translate_double_y();
      auto translate = shape_j.get_translate();
      auto rotation_j = shape_j.get_rotation();
      auto& nfp =
          comp_nfp(shape_j.base, rotation_j, shape_j.allowed_rotations,
                   shape_i.base, rotation_i, shape_i.allowed_rotations, layout);
      FT pd;
      // j < i时的pd已经被更新
      if (j > i) {
        Point_2 relative_point(ix - jx, iy - jy);

        pd = comp_pd(nfp, double_ix - double_jx, double_iy - double_jy, layout);
        layout.set_pd(i, j, pd);
      }
      c.nfps.push_back(&nfp);
      c.translations.push_back(translate);
      c.translate_x.push_back(double_jx);
    }
    c.initialize();
    auto points = c.get_perfect_points();
    // points为空时
    if (points.empty()) {
      std::cerr << "Error: Candidate points set is empty" << std::endl;
      break;
    }
    auto& first_point = points[0];
    auto first_x = first_point.x();
    auto first_y = first_point.y();
    auto double_first_x = CGAL::to_double(first_x);
    auto double_first_y = CGAL::to_double(first_y);
    for (size_t j = 0; j < num_poly; j++) {
      size_t k = j;
      if (i == j) {
        continue;
      } else if (j > i) {
        k--;
      }
      auto& nfp = c.nfps[k];
      // auto& t = c.translations[k];
      auto& shape_j = layout.sheet_parts[0][j];
      auto& jx = shape_j.get_translate_ft_x();
      auto& jy = shape_j.get_translate_ft_y();
      auto double_jx = shape_j.get_translate_double_x();
      auto double_jy = shape_j.get_translate_double_y();
      Point_2 relative_point(first_x - jx, first_y - jy);
      // 通常情况下此时pd为0
      auto pd = comp_pd(*nfp, double_first_x - double_jx,
                        double_first_y - double_jy, layout);
      shape_i.set_translate(first_x, first_y);
      layout.set_pd(i, j, pd);
    }
  }
  // 判断初始解是否可行
  auto pure_overlap = layout.get_pure_total_pd();
  if (pure_overlap > geo::BIAS) {
    throw std::runtime_error("Error get_init_solu(): initial solution is not feasible");
  }
  // 更新layout当前长度
  //s << "layout update" << std::endl;
  layout.update_cur_length();
  layout.best_length = layout.cur_length;
  // 更新sheet
  layout.sheets[0].set_width(layout.cur_length);
  //s << "update best result" << std::endl;
  layout.best_result = layout.sheet_parts;
  //s << "update best util" << std::endl;
  layout.best_utilization = CGAL::to_double(
      layout.area / (layout.best_length * layout.sheets[0].height));
  //s << "init solu end" << std::endl;
}
bool minimize_overlap(Layout& layout, volatile bool* requestQuit) {
  size_t numIterations = 0;
  FT minOverlap = geo::INF;
  layout.initialize_miu();
  std::clog << "It: ";
  std::vector<size_t> indices;
  for (size_t i = 0; i < layout.poly_num; i++) {
    indices.push_back(i);
  }
  while (numIterations < layout.maxIterations) {
    if (*requestQuit) {
      return false;
    }
    //  对polygon随机排列
    std::shuffle(indices.begin(), indices.end(), rand::random_engine3);
    // 随机选一个polygon，寻找它可能存在的最小位置
    for (size_t i = 0; i < layout.poly_num; i++) {
      auto idx = indices[i];
      auto& shape = layout.sheet_parts[0][idx];
      auto cur_pd = layout.get_one_polygon_pd(idx);
      // 如果这个polygon与其他polygon没有发生重叠，直接跳过
      if (cur_pd < geo::BIAS && cur_pd > -geo::BIAS) {
        continue;
      }
      // 遍历polygon的所有角度
      for (auto& rotation : shape.reduced_rotations) {
        auto rotate = geo::get_rotate(rotation, shape.allowed_rotations);
        Polygon_with_holes_2 rotated;
        if (rotate) {
          rotated = geo::transform_polygon_with_holes(*rotate, *shape.base);
        } else {
          rotated = *shape.base;
        }
        // 计算候选点
        CandidatePoints c(layout.poly_num - 1);
        c.set_boundary(geo::comp_ifr(layout.sheets[0].sheet, rotated));
        for (size_t k = 0; k < layout.poly_num; k++) {
          if (k == idx) {
            continue;
          }
          auto& shape_k = layout.sheet_parts[0][k];
          auto& nfp = comp_nfp(shape_k.base, shape_k.get_rotation(),
                               shape_k.allowed_rotations, shape.base, rotation,
                               shape.allowed_rotations, layout);
          auto translate_x = shape_k.get_translate_double_x();
          c.nfps.push_back(&nfp);
          c.translations.push_back(shape_k.get_translate());
          c.translate_x.push_back(translate_x);
        }
        c.initialize();
        auto points = c.get_arrangement_points();
        auto size = points.size();
        // std::clog << "The size of Candidate points is: " << size <<
        // std::endl;
        //  TODO 这里可以压缩到多线程求最小值的规约问题->求min(new_pd)
        for (size_t l = 0; l < size; ++l) {
          auto& point = points[l];
          auto point_x = point.x();
          auto point_y = point.y();
          auto double_point_x = CGAL::to_double(point_x);
          auto double_point_y = CGAL::to_double(point_y);
          double new_pd = 0;
          std::vector<FT> temp(layout.poly_num - 1);
          // 尝试将polygon idx平移到候选点上去，计算平移到哪个点polygon
          // idx与其他多边形的pd之和会降低
          for (size_t k = 0; k < layout.poly_num - 1; k++) {
            size_t m = k;
            if (k >= idx) {
              m = k + 1;
            }
            auto& nfp = c.nfps[k];
            auto& shape_m = layout.sheet_parts[0][m];
            auto double_mx = shape_m.get_translate_double_x();
            auto double_my = shape_m.get_translate_double_y();
            temp[k] = comp_pd(*nfp, double_point_x - double_mx,
                              double_point_y - double_my, layout);
            auto t = CGAL::to_double(temp[k]);
            new_pd += t * layout.get_miu(idx, m);
          }
          // 选择使得pd更小的那个点，如果不存在这样的点，则不会对layout发生任何变化
          if (new_pd < cur_pd) {
            shape.set(rotation, point_x, point_y);
            for (size_t k = 0; k < layout.poly_num - 1; k++) {
              if (k < idx) {
                layout.set_pd(idx, k, temp[k]);
              } else {
                layout.set_pd(idx, k + 1, temp[k]);
              }
            }
            cur_pd = new_pd;
          }
          // 如果这个候选点对其他polygon没有重叠，直接返回
          if (cur_pd < geo::BIAS) {
            break;
          }
        }
      }
    }
    auto pure_overlap = layout.get_pure_total_pd();
    // auto overlap = layout.get_total_pd();
    if (pure_overlap < geo::BIAS) {
      std::clog << std::endl;
      return true;
    } else if (pure_overlap < minOverlap) {
      std::clog << std::endl;
      std::clog << "Overlap: " << pure_overlap << std::endl;
      std::clog << "Total/Miss: " << layout.pd_count << "/" << layout.pd_miss
                << ", " << (double)layout.pd_miss / layout.pd_count
                << std::endl;
      std::clog << "It: ";
      minOverlap = pure_overlap;
      numIterations = 0;
    }
    layout.update_miu();
    numIterations++;
  }
  std::clog << std::endl;
  return false;
}
void GOMH(Layout& layout,
          size_t max_time,
          std::function<void(const Solution&)> ProgressHandler,
          volatile bool* requestQuit) {
  clock_t start = clock();
  get_init_solu(layout);
  ProgressHandler(Solution(
      CGAL::to_double(layout.best_length), layout.best_utilization,
      ((double)(clock() - start) / CLOCKS_PER_SEC), layout.best_result[0]));
  shrink(layout);
  while (((double)(clock() - start) / CLOCKS_PER_SEC) < max_time) {
    if (*requestQuit) {
      return;
    }
    auto feasible = minimize_overlap(layout, requestQuit);
    if (feasible) {
      layout.best_result = layout.sheet_parts;
      layout.best_length = (std::min)(layout.best_length, layout.cur_length);
      layout.best_utilization = CGAL::to_double(
          layout.area / (layout.best_length * layout.sheets[0].height));
      ProgressHandler(Solution(
          CGAL::to_double(layout.best_length), layout.best_utilization,
          ((double)(clock() - start) / CLOCKS_PER_SEC), layout.best_result[0]));
      // print_result(layout, start);
      if (layout.best_length <= layout.lower_length) {
        break;
      }
      shrink(layout);
    } else {
      auto t = layout.cur_length * (1 + layout.rinc);
      if (t >= layout.best_length) {
        shrink(layout);
      } else {
        expand(layout);
      }
    }
  }
  std::clog << "best length: " << layout.best_length << std::endl;
  std::clog << "best utilization rate: " << layout.best_utilization
            << std::endl;
}
}  // namespace nesting