bst_sub3.cxx

///////////////////////////////////////////////////////////////////////////////
// BSD 3-Clause License
//
// Copyright (c) 2018, The Regents of the University of California
// 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 copyright holder 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.
///////////////////////////////////////////////////////////////////////////////

/*
#**************************************************************************
#**   Please read README to see how to run this program
***   Created by Chung-Wen Albert Tsao  on May 2, 2000*
#**
#**
#**************************************************************************
*/

#include "bst_sub3.h"

#include "Global_var.h"
#include "bst.h"
#include "bst_header.h"
#include "bst_sub1.h"

#define NN 0
#define SS 1
#define EE 2
#define WW 3

int N_S = 0, N_S0 = 0;                 /* Number of detour points set  */
ObstacleType *Obstacle, *Tmp_Obstacle; /* which are rectangles */
double** ObsBox;
PointType *S, *S0; /* Set of detour points, which are obstacle vertices */
int* S_Par;
int* S_N_Hops;
int* In_Tree;
double* S_Cost;
TreeNodeType* TreeNode;
int BufferArea[2];

extern int on_line_segment(double, double, double, double, double, double);
extern int lineIntersect(PointType*,
                         PointType,
                         PointType,
                         PointType,
                         PointType);
extern int same_Point(PointType, PointType);
extern double Point_dist(PointType, PointType);
extern int
bbox_overlap(double, double, double, double, double, double, double, double);
extern int _bbox_overlap(double,
                         double,
                         double,
                         double,
                         double,
                         double,
                         double,
                         double,
                         double);
extern void build_trr(TrrType* ms, double d, TrrType* trr);
extern void make_intersect(TrrType* trr1, TrrType* trr2, TrrType* t);
extern void line2ms(TrrType* ms, PointType p1, PointType p2);
extern void core_mid_point(TrrType* trr, PointType* p);
extern int ms_type(TrrType* trr);
extern double ms_distance(TrrType* ms1, TrrType* ms2);
extern double pt2linedist(PointType p1,
                          PointType p2,
                          PointType p3,
                          PointType* ans);

extern double linedist(PointType lpt0,
                       PointType lpt1,
                       PointType lpt2,
                       PointType lpt3,
                       PointType pt[2]);
extern int Area_compare_inc(const void* p1, const void* q1);
extern int pt_on_line_segment(PointType pt, PointType pt1, PointType pt2);

extern int equal(double x, double y);
extern double pt_skew(PointType pt);
extern int calc_line_type(PointType pt1, PointType pt2);
extern void calc_vertices(AreaType* area);

extern void calc_pt_delays(AreaType* area,
                           PointType* q1,
                           PointType q0,
                           PointType q2);
extern void calc_BS_located(PointType* pt,
                            AreaType* area,
                            PointType* p1,
                            PointType* p2);
extern void calc_JS_delay(AreaType* area, AreaType* area_L, AreaType* area_R);

/********************************************************************/
/*                                                                  */
/********************************************************************/
double Point_E_dist(PointType p1, PointType p2)
{
  double a, b;

  a = ABS(p1.x - p1.x);
  b = ABS(p2.y - p2.y);
  return (sqrt(a * a + b * b));
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int cp_PointType_array(PointType a[], int n, PointType b[])
{
  int i;

  for (i = 0; i < n; ++i) {
    b[i] = a[i];
  }
  return (n);
}
/***********************************************************************/
/*    cp a[] to b[]                                                    */
/***********************************************************************/
int cp_area_array(AreaType a[], int n, AreaType b[])
{
  int i;

  for (i = 0; i < n; ++i) {
    b[i] = a[i];
  }
  return (n);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void append_area_array(AreaType a[], int n, AreaType b[], int* m)
{
  int i;

  for (i = 0; i < n; ++i) {
    b[*m + i] = a[i];
  }
  *m = *m + n;
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void check_area(AreaType* area)
{
  int i, j, k1, k2, n, line0type, line1type;
  PointType q0, q2, p1, p2;
  double a, b;

  assert(area->n_mr <= MAX_mr_PTS);
  if (Skew_B == 0) {
    assert(area->n_mr <= 2);
  }
  line0type = calc_line_type(area->line[0][0], area->line[0][1]);
  line1type = calc_line_type(area->line[1][0], area->line[1][1]);
  assert(line0type == line1type);

  n = area->n_mr;
  for (i = 0; i < n; ++i) {
    j = (i + 1) % n;
    q0 = area->mr[i];
    q2 = area->mr[j];
    a = ABS(q0.x - q2.x);
    b = ABS(q0.y - q2.y);
    if (!equal(a, 0) && !equal(b, 0)) {
      /* no skew reservation */
      if (equal(pt_skew(q0), Skew_B) == NO
          || equal(pt_skew(q2), Skew_B) == NO) {
        print_area(area);
        assert(0);
      }
    }
    if (line0type == MANHATTAN_ARC) {
      assert(equal(a, 0) || equal(b, 0) || equal(a, b));
    }
  }
  n = area->npts;
  for (i = 0; i < n; ++i) {
    k1 = area->vertex[i];
    k2 = area->vertex[(i + 1) % n];
    p1 = area->mr[k1];
    p2 = area->mr[k2];
    for (j = k1 + 1; j < k2; ++j) {
      assert(pt_on_line_segment(area->mr[j], p1, p2));
    }
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void check_area_array(AreaType area[], int n)
{
  int i;

  for (i = 0; i < n; ++i) {
    check_area(&(area[i]));
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int check_area_line(AreaType* area)
{
  int linetype;

  linetype = calc_line_type(area->line[0][0], area->line[0][1]);

  if (linetype == VERTICAL) {
    assert(area->line[0][0].y > area->line[0][1].y);
    assert(area->line[1][0].y > area->line[1][1].y);
  } else if (linetype == HORIZONTAL) {
    assert(area->line[0][0].x > area->line[0][1].x);
    assert(area->line[1][0].x > area->line[1][1].x);
  }
  return (linetype);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void check_all_areas_line(AreaType area[], int n)
{
  int i;

  for (i = 0; i < n; ++i) {
    check_area_line(&(area[i]));
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
double calc_subpath_length(PointType path[], int k1, int k2, int mode)
{
  int i;
  double cost;

  cost = 0;
  for (i = k1; i < k2; ++i) {
    if (mode == MANHATTAN) {
      cost += Point_dist(path[i], path[i + 1]);
    } else {
      cost += Point_E_dist(path[i], path[i + 1]);
    }
  }
  return (cost);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
double calc_pathlength(PointType path[], int n, int mode)
{
  double cost;

  cost = calc_subpath_length(path, 0, n - 1, mode);
  return (cost);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void reverse_mr(AreaType* area)
{
  int i, n;
  PointType tmp_Pt;

  n = area->n_mr;
  for (i = 0; i < n / 2; ++i) {
    tmp_Pt = area->mr[i];
    area->mr[i] = area->mr[n - 1 - i];
    area->mr[n - 1 - i] = tmp_Pt;
  }
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void get_block_bbox(double bbox[], ObstacleType* block)
{
  int i;

  bbox[EE] = bbox[NN] = -DBL_MAX;
  bbox[WW] = bbox[SS] = DBL_MAX;
  for (i = 0; i < block->n_vertex; ++i) {
    bbox[EE] = max(bbox[EE], block->vertex[i].x);
    bbox[WW] = min(bbox[WW], block->vertex[i].x);
    bbox[NN] = max(bbox[NN], block->vertex[i].y);
    bbox[SS] = min(bbox[SS], block->vertex[i].y);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_ms_JS(AreaType* area)
{
  int i, line0type, line1type;
  TrrType ms0, ms1, trr0, trr1, t0, t1;

  line0type = calc_line_type(area->line[0][0], area->line[0][1]);
  line1type = calc_line_type(area->line[1][0], area->line[1][1]);
  assert(line0type == line1type);
  assert(line0type == MANHATTAN_ARC);

  if (same_Point(area->line[0][0], area->line[0][1])
      && same_Point(area->line[1][0], area->line[1][1])) {
    return;
  }

  line2ms(&ms0, area->line[0][0], area->line[0][1]);
  line2ms(&ms1, area->line[1][0], area->line[1][1]);

  assert(equal(ms_distance(&ms0, &ms1), area->dist));

  build_trr(&ms0, area->dist, &trr0);
  build_trr(&ms1, area->dist, &trr1);
  make_intersect(&trr1, &ms0, &t0);
  make_intersect(&trr0, &ms1, &t1);
  ms2line(&t0, &(area->line[0][1]), &(area->line[0][0]));
  ms2line(&t1, &(area->line[1][1]), &(area->line[1][0]));
  assert(area->line[0][0].y >= area->line[0][1].y);
  assert(area->line[1][0].y >= area->line[1][1].y);
  assert(area->n_mr <= 2);
  if (area->n_mr == 2 && (area->L_EdgeLen == 0 || area->R_EdgeLen == 0)) {
    assert(equal(area->mr[0].max, area->mr[1].max));
    assert(equal(area->mr[0].min, area->mr[1].min));
    if (area->L_EdgeLen == 0) {
      i = 0;
    } else {
      i = 1;
    }
    assert(pt_on_line_segment(area->line[i][0], area->mr[0], area->mr[1]));
    assert(pt_on_line_segment(area->line[i][1], area->mr[0], area->mr[1]));
    area->mr[0].x = area->line[i][0].x;
    area->mr[0].y = area->line[i][0].y;
    if (same_Point(area->line[i][0], area->line[i][1])) {
      area->n_mr = area->npts = 1;
    } else {
      assert(area->npts == 2);
      area->mr[1].x = area->line[i][1].x;
      area->mr[1].y = area->line[i][1].y;
    }
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void get_JR_bbox(double bbox[], AreaType* area)
{
  int i, j;

  bbox[EE] = bbox[NN] = -DBL_MAX;
  bbox[WW] = bbox[SS] = DBL_MAX;

  for (i = 0; i < 2; ++i) {
    for (j = 0; j < 2; ++j) {
      bbox[EE] = max(bbox[EE], area->line[i][j].x);
      bbox[WW] = min(bbox[WW], area->line[i][j].x);
      bbox[NN] = max(bbox[NN], area->line[i][j].y);
      bbox[SS] = min(bbox[SS], area->line[i][j].y);
    }
  }
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void get_mr_bbox(double bbox[], AreaType* area)
{
  int i;

  bbox[EE] = bbox[NN] = -DBL_MAX;
  bbox[WW] = bbox[SS] = DBL_MAX;

  for (i = 0; i < area->n_mr; ++i) {
    bbox[EE] = max(bbox[EE], area->mr[i].x);
    bbox[WW] = min(bbox[WW], area->mr[i].x);
    bbox[NN] = max(bbox[NN], area->mr[i].y);
    bbox[SS] = min(bbox[SS], area->mr[i].y);
  }
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int JS_ms_type(AreaType* area)
{
  int i, j, line0type, line1type;
  TrrType ms0, ms1;
  PointType p0, p1;
  double x, y;

  line0type = calc_line_type(area->line[0][0], area->line[0][1]);
  line1type = calc_line_type(area->line[1][0], area->line[1][1]);

  assert(line0type == MANHATTAN_ARC && line1type == MANHATTAN_ARC);

  line2ms(&ms0, area->line[0][0], area->line[0][1]);
  line2ms(&ms1, area->line[1][0], area->line[1][1]);

  i = ms_type(&ms0);
  j = ms_type(&ms1);
  assert(i * j == 0 || i * j == 1);
  if (i == 0 && j == 0) {
    core_mid_point(&ms0, &p0);
    core_mid_point(&ms1, &p1);
    x = p0.x - p1.x;
    y = p0.y - p1.y;
    if (equal(x, 0) || equal(y, 0)) {
      return (0);
    } else if (x * y > 0) {
      return (-1);
    } else {
      return (1);
    }
  } else if (i != 0) {
    return (i);
  } else if (j != 0) {
    return (j);
  }
  assert(0);
  return (0);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void calc_N_S0()
{
  int i;

  N_S0 = 0;
  for (i = 0; i < N_Obstacle; ++i) {
    N_S0 += Obstacle[i].n_vertex;
  }

  S = (PointType*) calloc(N_S0 + 2, sizeof(PointType));
  assert(S != NULL);

  S0 = (PointType*) calloc(N_S0 + 2, sizeof(PointType));
  assert(S0 != NULL);

  S_Par = (int*) calloc(N_S0 + 2, sizeof(int));
  assert(S_Par != NULL);

  S_N_Hops = (int*) calloc(N_S0 + 2, sizeof(int));
  assert(S_N_Hops != NULL);

  In_Tree = (int*) calloc(N_S0 + 2, sizeof(int));
  assert(In_Tree != NULL);

  S_Cost = (double*) calloc(N_S0 + 2, sizeof(double));
  assert(S_Cost != NULL);
}

/***********************************************************************/
/* S[N_S] = set of obstacle vertices as detour points             */
/***********************************************************************/
void calc_S()
{
  int i, j, k, n;

  calc_N_S0();
  n = 0;
  for (i = 0; i < N_Obstacle; ++i) {
    k = Obstacle[i].n_vertex;
    for (j = 0; j < k; ++j) {
      S[n] = S0[n] = Obstacle[i].vertex[j];
      n++;
    }
  }
  assert(n == N_S0);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void init_S(PointType p1, PointType p2, double range)
{
  int i;
  double dist;

  N_S = 0;
  for (i = 0; i < N_S0; ++i) {
    dist = Point_dist(p1, S0[i]) + Point_dist(p2, S0[i]);
    if (dist <= range) {
      S[N_S++] = S0[i];
    }
  }
  S[N_S++] = p1;
  S[N_S++] = p2;
  assert(N_S <= N_S0 + 2);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void print_S()
{
  int i;
  FILE* f;

  f = fopen("t", "w");
  assert(f != NULL);

  fprintf(f, "\"S \n");
  for (i = 0; i < N_S; ++i) {
    fprintf(f, "move %f %f %d \n", S[i].x, S[i].y, i);
    fprintf(f, "     %f %f \n", S[i].x, S[i].y);
  }
  fprintf(f, "\n\"p1p2 \n");
  fprintf(f, "move %f %f \n", S[N_S - 2].x, S[N_S - 2].y);
  fprintf(f, "     %f %f \n", S[N_S - 1].x, S[N_S - 1].y);
  fclose(f);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int block_compare_inc1(const void* p1, const void* q1)
{
  int *p, *q;
  double a[4], b[4];

  p = (int*) p1;
  q = (int*) q1;

  get_block_bbox(a, &(Obstacle[*p]));
  get_block_bbox(b, &(Obstacle[*q]));
  return ((a[WW] > b[WW]) ? YES : NO);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int block_compare_inc2(const void* p1, const void* q1)
{
  int *p, *q;
  double a[4], b[4];

  p = (int*) p1;
  q = (int*) q1;

  get_block_bbox(a, &(Obstacle[*p]));
  get_block_bbox(b, &(Obstacle[*q]));
  return ((a[SS] > b[SS]) ? YES : NO);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void construct_TreeNode(int par,
                        int* index,
                        int obstacle[],
                        int n_obstacle,
                        int dir)
{
  int i, j, n1, n2, b[1000], cp;

  if (n_obstacle < 1) {
    return;
  }

  (*index)++;
  cp = *index;
  assert(cp >= 0 && cp < 2 * N_Obstacle - 1);
  TreeNode[cp].parent = par;

  if (par >= 0) {
    TreeNode[cp].level = TreeNode[par].level + 1;
    if (TreeNode[par].L == NIL) {
      TreeNode[par].L = cp;
    } else {
      assert(TreeNode[par].R == NIL);
      TreeNode[par].R = cp;
    }
  }
  TreeNode[cp].bbox[EE] = TreeNode[cp].bbox[NN] = -DBL_MAX;
  TreeNode[cp].bbox[WW] = TreeNode[cp].bbox[SS] = DBL_MAX;
  for (i = 0; i < n_obstacle; ++i) {
    /*
        get_block_bbox(obox, &(Obstacle[obstacle[i]]));
    */
    j = obstacle[i];
    TreeNode[cp].bbox[EE] = max(TreeNode[cp].bbox[EE], ObsBox[j][EE]);
    TreeNode[cp].bbox[NN] = max(TreeNode[cp].bbox[NN], ObsBox[j][NN]);
    TreeNode[cp].bbox[WW] = min(TreeNode[cp].bbox[WW], ObsBox[j][WW]);
    TreeNode[cp].bbox[SS] = min(TreeNode[cp].bbox[SS], ObsBox[j][SS]);
  }
  if (n_obstacle == 1) {
    TreeNode[cp].id = obstacle[0];
    if (CHECK) {
      printf("TreeNode[%d].level = %d \n", cp, TreeNode[cp].level);
    }
  } else {
    if (dir == 0) {
      qsort(obstacle, n_obstacle, sizeof(int), block_compare_inc1);
    } else {
      qsort(obstacle, n_obstacle, sizeof(int), block_compare_inc2);
    }

    n1 = n_obstacle / 2;
    n2 = n_obstacle - n1;
    assert(n2 < 1000);
    for (i = 0; i < n2; ++i) {
      b[i] = obstacle[n1 + i];
    }

    construct_TreeNode(cp, index, obstacle, n1, 1 - dir);
    construct_TreeNode(cp, index, b, n2, 1 - dir);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void init_construct_TreeNode()
{
  int i;
  int* a;

  i = -1;

  a = (int*) calloc(N_Obstacle, sizeof(int));
  assert(a != NULL);
  for (i = 0; i < N_Obstacle; ++i) {
    a[i] = i;
  }

  TreeNode = (TreeNodeType*) calloc(2 * N_Obstacle, sizeof(TreeNodeType));
  assert(TreeNode != NULL);

  for (i = 0; i < 2 * N_Obstacle; ++i) {
    TreeNode[i].L = TreeNode[i].R = TreeNode[i].parent = TreeNode[i].id = NIL;
    TreeNode[i].level = 0;
  }

  i = NIL;
  construct_TreeNode(NIL, &i, a, N_Obstacle, 0);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int unblocked_segment_sub1(PointType* p0, PointType* p1, int i)
{
  int j, k, ans;

  k = _bbox_overlap(p0->x,
                    p0->y,
                    p1->x,
                    p1->y,
                    TreeNode[i].bbox[WW],
                    TreeNode[i].bbox[SS],
                    TreeNode[i].bbox[EE],
                    TreeNode[i].bbox[NN],
                    0);

  if (k == NO) {
    return (YES);
  }

  j = TreeNode[i].id;
  if (j == NIL) { /* internal node */
    ans = unblocked_segment_sub1(p0, p1, TreeNode[i].L);
    if (ans == YES) {
      ans = unblocked_segment_sub1(p0, p1, TreeNode[i].R);
    }
  } else {
    assert(TreeNode[i].L == NIL && TreeNode[i].R == NIL);
    ans = 1 - line_into_rectangle(p0, p1, &(Obstacle[j]));
  }
  return (ans);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void read_obstacle_file_sub(FILE* f)
{
  int i, j, n;
  double x, y;

  fscanf(f, "%d \n", &N_Obstacle);

  Obstacle = (ObstacleType*) calloc(N_Obstacle, sizeof(ObstacleType));
  assert(Obstacle != NULL);

  Tmp_Obstacle = (ObstacleType*) calloc(N_Obstacle, sizeof(ObstacleType));
  assert(Tmp_Obstacle != NULL);

  ObsBox = (double**) calloc(N_Obstacle, sizeof(double*));
  assert(ObsBox != NULL);
  for (i = 0; i < N_Obstacle; i++) {
    ObsBox[i] = (double*) calloc(4, sizeof(double));
    assert(ObsBox[i] != NULL);
  }

  i = 0;
  while (fscanf(f, "%d:  \n", &n) == 1) {
    Obstacle[i].n_vertex = n;
    for (j = 0; j < n; ++j) {
      assert(2 == fscanf(f, "     %lf  %lf \n", &x, &y));
      Obstacle[i].vertex[j].x = x;
      Obstacle[i].vertex[j].y = y;
    }
    i++;
  }
  assert(i == N_Obstacle);

  for (i = 0; i < N_Obstacle; i++) {
    get_block_bbox(ObsBox[i], &(Obstacle[i]));
  }
  init_construct_TreeNode();
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void read_obstacle_file(const string& fn)
{
  FILE* f = fopen(fn.c_str(), "r");
  assert(f != NULL);
  read_obstacle_file_sub(f);
  calc_S();
  fclose(f);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int distinct_obstacle(int n)
{
  int i;
  double x1, x2, x3, x4, y1, y2, y3, y4;

  x1 = Obstacle[n].vertex[0].x;
  y1 = Obstacle[n].vertex[0].y;
  x2 = Obstacle[n].vertex[2].x;
  y2 = Obstacle[n].vertex[2].y;
  for (i = 0; i < n; ++i) {
    x3 = Obstacle[i].vertex[0].x;
    y3 = Obstacle[i].vertex[0].y;
    x4 = Obstacle[i].vertex[2].x;
    y4 = Obstacle[i].vertex[2].y;
    if (bbox_overlap(x1, y1, x2, y2, x3, y3, x4, y4)) {
      return (NO);
    }
  }
  return (YES);
}

/***************************************************************************/
/*  produce the the blockage.                                              */
/***************************************************************************/
void generate_obstacle_file(char fn[], int n_blocks, int max_coor)
{
  FILE* out;
  double x, y, w, h;
  int i, j, max_block_size;
  char a[100], *b = ".block";

  N_Obstacle = n_blocks;
  Obstacle = (ObstacleType*) calloc(N_Obstacle, sizeof(ObstacleType));
  assert(Obstacle != NULL);

  strcpy(a, fn);
  strcat(a, b);

  out = fopen(a, "w");
  assert(out != NULL);

  fprintf(out, "%d  \n", n_blocks);
  srand48(1);
  max_block_size = max_coor / 20;
  i = 0;
  while (i < n_blocks) {
    x = rand() % max_coor;
    y = rand() % max_coor;
    w = max_block_size / 6 + rand() % max_block_size;
    h = max_block_size / 6 + rand() % max_block_size;
    Obstacle[i].n_vertex = 4;
    Obstacle[i].vertex[0].x = x - w;
    Obstacle[i].vertex[1].x = x - w;
    Obstacle[i].vertex[2].x = x + w;
    Obstacle[i].vertex[3].x = x + w;
    Obstacle[i].vertex[0].y = y - h;
    Obstacle[i].vertex[1].y = y + h;
    Obstacle[i].vertex[2].y = y + h;
    Obstacle[i].vertex[3].y = y - h;
    if (distinct_obstacle(i)) {
      fprintf(out, "4: \n");
      for (j = 0; j < 4; ++j) {
        x = Obstacle[i].vertex[j].x;
        y = Obstacle[i].vertex[j].y;
        fprintf(out, "    %f  %f \n", x, y);
      }
      ++i;
    }
  }
  fclose(out);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int pt_outside_a_rectangle(PointType* pt, ObstacleType* block)
{
  double bbox[4];

  get_block_bbox(bbox, block);
  if (equal(pt->x, bbox[EE]) || equal(pt->x, bbox[WW])) {
    if (pt->y <= bbox[NN] + FUZZ && pt->y >= bbox[SS] - FUZZ) {
      return (ON);
    } else {
      return (OUTSIDE);
    }
  }
  if (equal(pt->y, bbox[NN]) || equal(pt->y, bbox[SS])) {
    if (pt->x <= bbox[EE] + FUZZ && pt->x >= bbox[WW] - FUZZ) {
      return (ON);
    } else {
      return (OUTSIDE);
    }
  }
  if (pt->x <= bbox[EE] && pt->x >= bbox[WW] && pt->y <= bbox[NN]
      && pt->y >= bbox[SS]) {
    return (INSIDE);
  } else {
    return (OUTSIDE);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int pt_outside_rectangles(PointType* pt)
{
  int i, j;

  for (i = 0; i < N_Obstacle; ++i) {
    j = pt_outside_a_rectangle(pt, &(Obstacle[i]));
    if (j == 0 || j == -1) {
      return (NO);
    }
  }
  return (YES);
}

/***************************************************************************/
/*  produce the test cases.                                                */
/***************************************************************************/
static void generate_terminals(char fn[], int max_coor, int n_sinks)
{
  FILE* out;
  int i;
  PointType pt;

  generate_obstacle_file(fn, 40, max_coor);

  out = fopen(fn, "w");
  assert(out != NULL);

  long current_time;
  time(&current_time);
  // string timestamp = U_Printtf("%s",ctime(&current_time));
  string timestamp = ctime(&current_time);

  string msg = "# UCLA clock benchmark 1.0";
  msg += "# Created       : " + timestamp;
  msg += "\n";
  msg += "# User          : tsao@cs.ucla.edu";
  msg += "\n";
  msg += "# PlatForm      : SunOS 5.7 sparc SUNW,Ultra-1";
  msg += "\n";
  msg += "# Source        : Randomly generated testcases ";
  msg += "\n";
  msg += "# Note          : Coordinate unit can be micro-meter or anything";
  msg += "\n";

  fprintf(out, "%s  \n\n", msg.c_str());
  fprintf(out, "NumPins : %d  \n\n", n_sinks);
  fprintf(out, "PerUnitResistance : 0.016600  Ohm \n\n");
  fprintf(out, "PerUnitCapacitance : .27000e-17  Farad \n\n");
  // fprintf(out,"%d  0.016600  2.7000e-17 \n", n_sinks);
  srand48(1);
  i = 0;
  while (i < n_sinks) {
    pt.x = rand() % max_coor;
    pt.y = rand() % max_coor;
    if (pt_outside_rectangles(&pt)) {
      fprintf(out, "Sink : %d  \n", i);
      fprintf(out, "     Coordinate : %d  %d \n", (int) pt.x, (int) pt.y);
      fprintf(out, "     Capacitive Load  : 5.0e-13\n");
      fprintf(out, "     Downstream_Delay : 0.0e-12\n");
      // fprintf(out,"%d  %d  %d 5.0e-13\n",i, (int) pt.x,(int) pt.y);
      i++;
    }
  }
  fclose(out);
  exit(0);
}

/**********************************************************************/
/* Copyright (C) 1994-2000  by Andrew B. Kahng, C.-W. Albert Tsao    */
/**********************************************************************/
bool BST_GenTestcase::GenerateTestcase(unsigned i)
{
  if (i > 0) {
    char fn[100];
    sprintf(fn, "n%d", i);
    generate_terminals(fn, 1000, i);
  }

  return i > 0;
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int distinct_Point(PointType q[], int n)
{
  int i;

  for (i = 0; i < n; ++i) {
    if (same_Point(q[n], q[i])) {
      return (NO);
    }
  }
  return (YES);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int line_x_rectangle(PointType* p0,
                     PointType* p1,
                     PointType q[],
                     ObstacleType* block)
{
  PointType v0, v1;
  int ans;
  int i, n, n_xing = 0;

  n = block->n_vertex;
  for (i = 0; i < n; ++i) {
    v0 = block->vertex[i];
    v1 = block->vertex[(i + 1) % n];
    ans = lineIntersect(&(q[n_xing]), *p0, *p1, v0, v1);
    if (ans == SAME_EDGE || ans == OVERLAPPING) {
      return (0);
    } else if (ans > 0 && distinct_Point(q, n_xing)) {
      q[n_xing].max = p0->max;
      q[n_xing].min = p0->min;
      n_xing++;
    }
  }
  assert(n_xing <= 2);
  return (n_xing);
}
/***********************************************************************/
/*  check if line p0,p1 goes into rectangle[k]                         */
/***********************************************************************/
int line_into_rectangle(PointType* p0, PointType* p1, ObstacleType* block)
{
  int ans, loc0, loc1, n_xing;
  PointType q[10];

  ans = bbox_overlap(p0->x,
                     p0->y,
                     p1->x,
                     p1->y,
                     block->vertex[0].x,
                     block->vertex[0].y,
                     block->vertex[2].x,
                     block->vertex[2].y);

  if (ans == NO) {
    return (NO);
  }
  loc0 = pt_outside_a_rectangle(p0, block);
  loc1 = pt_outside_a_rectangle(p1, block);

  if (loc0 == -1 || loc1 == -1) { /* p0 or p1 is inside rectangle */
    return (YES);
  } else if (same_Point(*p0, *p1)) { /* p0, p1 either on boundary or outside */
    return (NO);
  } else {
    n_xing = line_x_rectangle(p0, p1, q, block);
    if (n_xing == 2) {
      return (YES);
    } else {
      return (NO);
    }
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int unblocked_segment_sub2(PointType* p0, PointType* p1)
{
  int i;
  double bbox[4];

  for (i = 0; i < N_Obstacle; ++i) {
    get_block_bbox(bbox, &(Obstacle[i]));
    if (_bbox_overlap(p0->x,
                      p0->y,
                      p1->x,
                      p1->y,
                      bbox[WW],
                      bbox[SS],
                      bbox[EE],
                      bbox[NN],
                      FUZZ)) {
      if (line_into_rectangle(p0, p1, &(Obstacle[i]))) {
        return (NO);
      }
    }
  }
  return (YES);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int unblocked_segment(PointType* p0, PointType* p1)
{
  int ans;

  if (N_Obstacle == 0) {
    return (YES);
  }

  /*
    ans = unblocked_segment_sub2(p0,p1);
  */
  ans = unblocked_segment_sub1(p0, p1, 0);
  return (ans);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int same_test_pt(PointType p1, PointType p2, PointType p[3][2])
{
  int i;

  for (i = 0; i < 3; ++i) {
    if (same_Point(p1, p[i][0]) && same_Point(p2, p[i][1])) {
      return (YES);
    }
    if (same_Point(p1, p[i][1]) && same_Point(p2, p[i][0])) {
      return (YES);
    }
  }
  return (NO);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int legal_detour_path(PointType p1, PointType p2)
{
  int ans;
  PointType p[3][2];

  p[0][0].x = 18;
  p[0][0].y = 334;
  p[0][1].x = 408;
  p[0][1].y = 275;
  p[1][0] = p[0][1];
  p[1][1].x = 436;
  p[1][1].y = 275;
  p[2][0] = p[1][1];
  p[2][1].x = 476;
  p[2][1].y = 344;
  if (same_Point(p1, p2)) {
    return (NO);
  }

  ans = unblocked_segment(&p1, &p2);
  if (same_test_pt(p1, p2, p)) {
    assert(ans == YES);
  }

  return (ans);
}

/***********************************************************************/
/*  update the distance of a node i outside SPT to the new tree node cp*/
/***********************************************************************/
void update_free_node_info(int cp, int i)
{
  double cost;

  assert(In_Tree[cp] == YES); /* S[cp] is newly added into SPT */
  assert(In_Tree[i] == NO);   /* S[i] is  outside SPT */

  cost = S_Cost[cp] + Point_dist(S[cp], S[i]);
  if ((equal(cost, S_Cost[i]) && S_N_Hops[i] > S_N_Hops[cp] + 1)
      || (cost < S_Cost[i])) {
    S_Par[i] = cp;
    S_Cost[i] = cost;
    S_N_Hops[i] = S_N_Hops[cp] + 1;
  }
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void init_path_finder(int source, int dest)
{
  int i;

  for (i = 0; i < N_S0 + 2; ++i) {
    In_Tree[i] = NO;
    S_Par[i] = NIL;
    S_Cost[i] = DBL_MAX;
    S_N_Hops[i] = 0;
  }
  In_Tree[source] = YES;
  S_Par[source] = NIL;
  S_Cost[source] = 0;

  assert(same_Point(S[source], S[dest]) == NO);
  assert(legal_detour_path(S[source], S[dest]) == NO);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void print_path(FILE* f, PointType path[], int n)
{
  int i;

  fprintf(f, "move %.9f %.9f  n_hop:%d\n", path[0].x, path[0].y, n);
  for (i = 1; i < n; ++i) {
    fprintf(f, "draw %.9f %.9f  \n", path[i].x, path[i].y);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void shrink_an_obstacle(ObstacleType* block, double s)
{
  int i, n;
  double cx, cy;

  n = block->n_vertex;
  cx = cy = 0;
  for (i = 0; i < n; ++i) {
    cx += block->vertex[i].x;
    cy += block->vertex[i].y;
  }
  cx /= n;
  cy /= n;
  for (i = 0; i < n; ++i) {
    block->vertex[i].x = (cx + (s - 1) * block->vertex[i].x) / s;
    block->vertex[i].y = (cy + (s - 1) * block->vertex[i].y) / s;
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void print_an_obstacle(FILE* f, ObstacleType block, int index)
{
  int i, n;
  PointType pt;
  double bbox[4];

  if (1) {
    shrink_an_obstacle(&block, 7);
  }
  get_block_bbox(bbox, &block);

  n = block.n_vertex;

  pt = block.vertex[0];
  fprintf(f, "move %.9f %.9f  %d ", pt.x, pt.y, index);
  if (bbox[WW] > 260 && bbox[EE] < 350 && bbox[NN] < 240 && bbox[SS] > 100) {
    printf("   here");
  }
  fprintf(f, "\n");
  for (i = 0; i < n; ++i) {
    pt = block.vertex[(i + 1) % n];
    fprintf(f, "draw %.9f %.9f  \n", pt.x, pt.y);
  }
  for (i = 0; i < n / 2; ++i) {
    pt = block.vertex[i];
    fprintf(f, "move %.9f %.9f  \n", pt.x, pt.y);
    pt = block.vertex[(i + n / 2) % n];
    fprintf(f, "draw %.9f %.9f  \n", pt.x, pt.y);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void print_obstacles_sub(FILE* f)
{
  int i;

  for (i = 0; i < N_Obstacle; ++i) {
    print_an_obstacle(f, Obstacle[i], i);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void print_obstacles(const string& fn)
{
  string xfn = fn + ".xg";

  FILE* f = fopen(xfn.c_str(), "w");
  assert(f != NULL);
  print_obstacles_sub(f);
  fclose(f);
}

/***********************************************************************/
/*  calc the shortest path[] from the shortest path tree in v[]        */
/***********************************************************************/
int get_path_from_SPT(PointType path[])
{
  int i, n;
  PointType tmpPt;

  n = 0;
  i = N_S - 1;           /* S[N_S-1] is the destination */
  while (i != N_S - 2) { /* S[NS-2] is the source */
    path[n++] = S[i];
    i = S_Par[i];
  }
  path[n++] = S[i];
  for (i = 0; i < n / 2; ++i) {
    tmpPt = path[i];
    path[i] = path[n - 1 - i];
    path[n - 1 - i] = tmpPt;
  }

  return (n);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int path_finder_sub2(int cp)
{
  int i, min_i;
  double min_cost;

  In_Tree[cp] = YES; /* add S[cp] into SPT */
  for (i = 0; i < N_S; ++i) {
    if (In_Tree[i] == NO && legal_detour_path(S[cp], S[i])) {
      update_free_node_info(cp, i);
    }
  }
  min_cost = DBL_MAX;
  min_i = NIL; /* S[min_i] is next node to be added into SPT */
  for (i = 0; i < N_S; ++i) {
    if (In_Tree[i] == NO && S_Cost[i] < DBL_MAX) {
      if (equal(min_cost, S_Cost[i])) {
        if (S_N_Hops[i] < S_N_Hops[min_i]) {
          min_cost = S_Cost[i];
          min_i = i;
        }
      } else if (S_Cost[i] < min_cost) {
        min_cost = S_Cost[i];
        min_i = i;
      }
    }
  }
  if (N_S >= N_S0 + 2) {
    assert(min_i != NIL);
  }
  return (min_i); /* expand the node in S with minimum cost */
}

/***********************************************************************/
/*   find a planar shortest path connecting source and dest            */
/***********************************************************************/
double path_finder_sub(int source, int dest, PointType path[], int* n_hops)
{
  int cp;
  double cost;

  if (0 == 1) {
    printf("*************************************************\n");
    printf("source = S[%d](%f, %f), dest = S[%d](%f, %f) \n",
           source,
           S[source].x,
           S[source].y,
           dest,
           S[dest].x,
           S[dest].y);
    printf("*************************************************\n");
  }
  *n_hops = 0;
  init_path_finder(source, dest);
  cp = source;
  while (cp != dest) { /* S[dest] not in Shortest Path Tree yet */
    cp = path_finder_sub2(cp);
    if (cp == NIL) {
      return (DBL_MAX);
    }
  }
  In_Tree[cp] = YES;
  *n_hops = get_path_from_SPT(path);
  cost = calc_pathlength(path, *n_hops, MANHATTAN);
  assert(equal(S_Cost[dest], cost));
  return (cost);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
double path_finder(PointType p1, PointType p2, PointType path[], int* n_hops)
{
  double cost, range, dist, t;

  range = cost = dist = Point_dist(p1, p2);

  if (N_Obstacle == 0 || same_Point(p1, p2) || legal_detour_path(p1, p2)) {
    path[0] = p1;
    path[1] = p2;
    *n_hops = 2;
  } else {
    cost = DBL_MAX;
    do {
      if (cost == DBL_MAX) {
        range = range + dist;
      } else {
        range = cost;
      }
      init_S(p1, p2, range);
      t = path_finder_sub(N_S - 2, N_S - 1, path, n_hops);
      cost = t;
      if (N_S >= N_S0 + 2 && *n_hops <= 2) {
        print_S();
        assert(0);
      }
    } while (cost > range + FUZZ);
  }

  assert(*n_hops >= 2);
  assert(*n_hops <= 100);
  assert(cost <= INT_MAX);
  return (cost);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void test_path_finder_sub(FILE* f, PointType p0, PointType p1)
{
  PointType path[100];
  int n;
  double cost;

  cost = path_finder(p0, p1, path, &n);
  fprintf(f, "\"path(%.0f) \n", cost);
  print_path(f, path, n);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void test_path_finder(FILE* f)
{
  PointType p0, p1;

  p0.x = 448;
  p0.y = 334;
  p1.x = 408;
  p1.y = 343;
  test_path_finder_sub(f, p0, p1);
  p0.x = 408;
  p0.y = 275;
  p1.x = 436;
  p1.y = 275;
  test_path_finder_sub(f, p0, p1);
  p0.x = 436;
  p0.y = 275;
  p1.x = 476;
  p1.y = 344;
  test_path_finder_sub(f, p0, p1);
  p0.x = 18;
  p0.y = 334;
  p1.x = 476;
  p1.y = 344;
  test_path_finder_sub(f, p0, p1);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void calc_n_vertex(AreaType* area)
{
  area->vertex[0] = 0;
  area->vertex[1] = 1;
  if (same_Point(area->mr[0], area->mr[1])) {
    area->npts = area->n_mr = 1;
  } else {
    area->npts = area->n_mr = 2;
  }
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void get_mr_for_line_sub(AreaType* area, AreaType* new_area, int begin, int end)
{
  int i, k;

  k = new_area->n_mr;
  for (i = begin; i <= end; ++i) {
    new_area->mr[k++] = area->mr[i];
  }
  new_area->n_mr = k;
  new_area->vertex[0] = 0;

  for (i = 0; i < area->n_mr - 1; ++i) {
    assert(!same_Point(area->mr[i], area->mr[i + 1]));
  }

  if (new_area->n_mr > 1) {
    new_area->npts = 2;
    new_area->vertex[1] = new_area->n_mr - 1;
  } else {
    new_area->npts = 1;
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void get_mr_for_line(PointType pt, AreaType* area, AreaType* new_area, int dir)
{
  int i;

  for (i = 0; i < area->n_mr; ++i) {
    if (same_Point(pt, area->mr[i])) {
      new_area->n_mr = 0;
      if (dir == 0) {
        get_mr_for_line_sub(area, new_area, 0, i);
      } else {
        get_mr_for_line_sub(area, new_area, i, area->n_mr - 1);
      }
      assert(new_area->n_mr <= 4);
      return;
    }
  }
  for (i = 0; i < area->n_mr - 1; ++i) {
    if (pt_on_line_segment(pt, area->mr[i], area->mr[i + 1])) {
      calc_pt_delays(area, &pt, area->mr[i], area->mr[i + 1]);
      if (dir == 0) {
        new_area->n_mr = 0;
        get_mr_for_line_sub(area, new_area, 0, i);
        new_area->mr[new_area->n_mr] = pt;
        new_area->vertex[1] = new_area->n_mr;
        new_area->n_mr += 1;
        new_area->npts = 2;
      } else {
        new_area->mr[0] = pt;
        new_area->n_mr = 1;
        get_mr_for_line_sub(area, new_area, i + 1, area->n_mr - 1);
      }
      assert(new_area->n_mr <= 4);
      return;
    }
  }
  assert(0 == 1);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int modify_area_sub3(AreaType* area,
                     AreaType tmpArea[],
                     int* n,
                     ObstacleType* block)
{
  PointType p0, p1, q[10];
  int n_xing = 0;

  p0 = area->mr[0];
  p1 = area->mr[area->n_mr - 1];
  /* p0 and p1 are either outside or on the boundary of obstacle */

  n_xing = line_x_rectangle(&p0, &p1, q, block);
  if (n_xing == 2) {
    if (Point_dist(p0, q[0]) < Point_dist(p0, q[1])) {
      get_mr_for_line(q[0], area, &(tmpArea[*n]), 0);
      get_mr_for_line(q[1], area, &(tmpArea[*n + 1]), 1);
    } else {
      get_mr_for_line(q[1], area, &(tmpArea[*n]), 0);
      get_mr_for_line(q[0], area, &(tmpArea[*n + 1]), 1);
    }
    *n += 2;
    if (0) {
      check_area_array(tmpArea, *n);
    }
  } else {
    tmpArea[(*n)++] = *area;
  }
  return (n_xing);
}
/***********************************************************************/
/*  check if line p0,p1 goes into rectangle[k]                         */
/***********************************************************************/
void modify_area_sub2(AreaType* area,
                      AreaType tmpArea[],
                      int* n,
                      ObstacleType* block)
{
  int i, n_mr, loc0, loc1, n_xing;
  PointType p0, p1, q[10];

  n_mr = area->n_mr;
  loc0 = pt_outside_a_rectangle(&(area->mr[0]), block);
  loc1 = pt_outside_a_rectangle(&(area->mr[n_mr - 1]), block);
  if (loc0 > loc1) {
    reverse_mr(area);
    i = loc0;
    loc0 = loc1;
    loc1 = i;
  }
  p0 = area->mr[0];
  p1 = area->mr[n_mr - 1];
  assert(!same_Point(p0, p1));

  if (loc0 == INSIDE) {
    if (loc1 == ON) {
      tmpArea[*n].npts = tmpArea[*n].n_mr = 1;
      tmpArea[*n].mr[0] = p1;
      (*n)++;
    } else if (loc1 == OUTSIDE) {
      n_xing = line_x_rectangle(&p0, &p1, q, block);
      assert(n_xing == 1);
      get_mr_for_line(q[0], area, &(tmpArea[*n]), 1);
      (*n)++;
    }
  } else {
    /* p0 and p1 are either outside or on the boundary of obstacle */
    modify_area_sub3(area, tmpArea, n, block);
  }
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_area_sub1(AreaType* area,
                      AreaType tmpArea[],
                      int* n,
                      ObstacleType* block)
{
  int type;

  if (0) {
    check_area(area);
  }
  if (area->n_mr == 1) {
    type = pt_outside_a_rectangle(&(area->mr[0]), block);
    if (type == 1 || type == 0) { /* add one area */
      tmpArea[(*n)++] = *area;
    }
  } else { /* area->n_mr==2 */
    tmpArea[*n] = tmpArea[*n + 1] = *area;
    modify_area_sub2(area, tmpArea, n, block);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void expand_a_dimention_sub(int index, double z, int dir)
{
  double bbox[4];

  assert(index < N_Obstacle);
  get_block_bbox(bbox, &(Tmp_Obstacle[index]));
  if (dir == NN || dir == EE) {
    bbox[dir] = max(z + 1E8, bbox[dir]);
  } else {
    bbox[dir] = min(z - 1E8, bbox[dir]);
  }
  Tmp_Obstacle[index].vertex[0].x = bbox[WW];
  Tmp_Obstacle[index].vertex[1].x = bbox[WW];
  Tmp_Obstacle[index].vertex[2].x = bbox[EE];
  Tmp_Obstacle[index].vertex[3].x = bbox[EE];
  Tmp_Obstacle[index].vertex[0].y = bbox[SS];
  Tmp_Obstacle[index].vertex[1].y = bbox[NN];
  Tmp_Obstacle[index].vertex[2].y = bbox[NN];
  Tmp_Obstacle[index].vertex[3].y = bbox[SS];
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int bblock_overlap(ObstacleType* b1, ObstacleType* b2)
{
  double bbox1[4], bbox2[4];

  get_block_bbox(bbox1, b1);
  get_block_bbox(bbox2, b2);

  if (bbox1[SS] >= bbox2[NN]) {
    return (NO);
  }
  if (bbox2[SS] >= bbox1[NN]) {
    return (NO);
  }
  if (bbox1[WW] >= bbox2[EE]) {
    return (NO);
  }
  if (bbox2[WW] >= bbox1[EE]) {
    return (NO);
  }
  return (YES);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int expand_a_dimention_sub2(double z,
                            int dir,
                            int* n_unexp,
                            int unexp[],
                            int expanded[])
{
  int i, j, x, y, change = NO, n_expanded, n;

  n = *n_unexp;
  n_expanded = N_Obstacle - n;
  for (i = 0; i < n; ++i) {
    x = unexp[i];
    for (j = 0; j < n_expanded; ++j) {
      y = expanded[j];
      if (bblock_overlap(&(Tmp_Obstacle[x]), &(Tmp_Obstacle[y]))) {
        change = YES;
        expand_a_dimention_sub(x, z, dir);
        unexp[i] = unexp[*n_unexp - 1];
        (*n_unexp)--;
        expanded[n_expanded++] = x;
        return (YES);
      }
    }
  }
  return (change);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void expand_a_dimention(int index, double z, int dir)
{
  int i, unexp[1000], expanded[1000];
  int n_unexp, change;

  expand_a_dimention_sub(index, z, dir);

  for (i = 0; i < N_Obstacle; ++i) {
    unexp[i] = i;
  }
  unexp[index] = N_Obstacle - 1;
  n_unexp = N_Obstacle - 1;
  expanded[0] = index;

  do {
    change = expand_a_dimention_sub2(z, dir, &n_unexp, unexp, expanded);
  } while (change == YES);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int expand_obstacle(AreaType* area, int index, double bbox[4], int type)
{
  int i, j, k, ans, D0, D1;
  PointType p0, p1;

  if (type == 0) {
    if (same_Point(area->line[0][0], area->line[1][0])) {
      return (NO);
    } else {
      ans = line_into_rectangle(
          &(area->line[0][0]), &(area->line[1][0]), &(Tmp_Obstacle[index]));
      return (ans);
    }
  } else {
    if (type == 1) {
      D0 = EE;
      D1 = WW;
    } else {
      D0 = WW;
      D1 = EE;
    }
    p0.x = bbox[D0];
    p0.y = bbox[NN];
    p1.x = bbox[D1];
    p1.y = bbox[SS];
    i = line_into_rectangle(&p0, &(area->line[0][0]), &(Tmp_Obstacle[index]));
    j = line_into_rectangle(&p0, &(area->line[1][0]), &(Tmp_Obstacle[index]));
    if (i || j) {
      expand_a_dimention(index, bbox[NN], NN);
      expand_a_dimention(index, bbox[D0], D0);
      k = 1;
    } else {
      k = 0;
    }
    i = line_into_rectangle(&p1, &(area->line[0][1]), &(Tmp_Obstacle[index]));
    j = line_into_rectangle(&p1, &(area->line[1][1]), &(Tmp_Obstacle[index]));
    if (i || j) {
      if (k == 1) {
        return (YES);
      }
      expand_a_dimention(index, bbox[SS], SS);
      expand_a_dimention(index, bbox[D1], D1);
    }
    return (NO);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_area_sub0(AreaType area_out[], int* n_out)
{
  int n_in, i, j;
  AreaType area_in[10], area;

  area = area_out[0];
  *n_out = 1;
  for (i = 0; i < N_Obstacle; ++i) {
    n_in = cp_area_array(area_out, *n_out, area_in);
    *n_out = 0;
    for (j = 0; j < n_in; ++j) {
      modify_area_sub1(&(area_in[j]), area_out, n_out, &(Tmp_Obstacle[i]));
    }
    if (*n_out == 0) {
      return;
    }
  }
  assert(*n_out < 10);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
double path_between_JSline_sub(AreaType* area,
                               PointType line[2][2],
                               PointType path[100],
                               int* n,
                               PointType pt[])
{
  int done = NO, count = 0, min_n;
  PointType q0, q1, min_path[100];
  double cost, min_cost = DBL_MAX;

  assert(pt_on_line_segment(pt[0], line[0][0], line[0][1]));
  assert(pt_on_line_segment(pt[1], line[1][0], line[1][1]));

  done = NO;
  min_cost = DBL_MAX;
  while (!done) {
    count++;
    if (area != NULL) {
      assert(area->n_mr <= MAX_mr_PTS);
    }
    cost = path_finder(pt[0], pt[1], path, n);
    if (area != NULL) {
      assert(area->n_mr <= MAX_mr_PTS);
    }
    assert(*n <= 100);
    pt2linedist(path[1], line[0][0], line[0][1], &q0);
    pt2linedist(path[*n - 2], line[1][0], line[1][1], &q1);
    if (!same_Point(pt[0], q0) || !same_Point(pt[1], q1)) {
      pt[0] = q0;
      pt[1] = q1;
    } else {
      done = YES;
    }
    if (cost < min_cost) {
      min_cost = cost;
      min_n = cp_PointType_array(path, *n, min_path);
    } else {
      done = YES;
    }
  }
  if (area != NULL) {
    assert(area->n_mr <= MAX_mr_PTS);
  }
  *n = cp_PointType_array(min_path, min_n, path);
  assert(pt_on_line_segment(path[0], line[0][0], line[0][1]));
  assert(pt_on_line_segment(path[*n - 1], line[1][0], line[1][1]));
  return (min_cost);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void center_of_segment(PointType* q, PointType p1, PointType p2)
{
  q->x = (p1.x + p2.x) / 2;
  q->y = (p1.y + p2.y) / 2;
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
double path_between_JSline(AreaType* area,
                           PointType line[2][2],
                           PointType path[100],
                           int* n)
{
  int i, j, m, num[2], min_n;
  PointType q[2][9], pt[2], min_path[100];
  double cost, min_cost;

  linedist(line[0][0], line[0][1], line[1][0], line[1][1], pt);

  for (i = 0; i < 2; ++i) {
    m = 0;
    q[i][m++] = pt[i];
    if (!same_Point(line[i][0], line[i][1])) {
      q[i][m++] = line[i][0];
      q[i][m++] = line[i][1];
      center_of_segment(&(q[i][m++]), line[i][0], line[i][1]);
    }
    num[i] = m;
  }

  num[0] = num[1] = 1;
  min_cost = DBL_MAX;
  for (i = 0; i < num[0]; ++i) {
    pt[0] = q[0][i];
    for (j = 0; j < num[1]; ++j) {
      pt[1] = q[1][j];
      cost = path_between_JSline_sub(area, line, path, n, pt);
      assert(*n <= 100);
      if (cost < min_cost) {
        min_cost = cost;
        min_n = cp_PointType_array(path, *n, min_path);
      }
    }
  }
  *n = cp_PointType_array(min_path, min_n, path);
  return (min_cost);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int get_boundary_segment(AreaType* area, AreaType new_area[])
{
  int j, k, m, n, npts;

  n = 0;
  if (0) {
    calc_vertices(area);
  }
  npts = area->npts;
  if (npts > 2) {
    for (j = 0; j < npts; ++j) {
      new_area[n] = *area;
      if (0) {
        check_area(&(new_area[n]));
      }
      new_area[n].npts = 2;
      m = 0;
      k = area->vertex[j] - 1;
      do {
        k = (k + 1) % area->n_mr;
        new_area[n].mr[m++] = area->mr[k];
      } while (k != area->vertex[(j + 1) % npts]);
      new_area[n].n_mr = m;
      new_area[n].vertex[0] = 0;
      new_area[n].vertex[1] = m - 1;
      calc_vertices(&(new_area[n]));
      assert(new_area[n].npts == 2);
      if (0) {
        check_area(&(new_area[n]));
      }
      n++;
    }
    assert(m <= 4);
  } else {
    new_area[n] = *area;
    modify_ms_JS(&(new_area[n]));
    n++;
  }
  return (n);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int check_lines(AreaType* area)
{
  assert(area->line[0][0].y >= area->line[0][1].y);
  assert(area->line[1][0].y >= area->line[1][1].y);
  if (same_Point(area->line[0][0], area->line[1][0])
      && same_Point(area->line[0][1], area->line[1][1])) {
    return (NO);
  } else {
    return (YES);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int modify_area_case1(AreaType* area, AreaType out_area[])
{
  int i, tmp_n, out_n, in_n, containing, type;
  AreaType in_area[6], tmp_area[100];
  double bbox[4];

  in_n = get_boundary_segment(area, in_area);
  assert(in_n <= 6);

  get_JR_bbox(bbox, &(in_area[0]));
  type = JS_ms_type(&(in_area[0]));
  check_lines(&(in_area[0]));

  for (i = 0; i < N_Obstacle; ++i) {
    Tmp_Obstacle[i] = Obstacle[i];
  }
  for (i = 0; i < N_Obstacle; ++i) {
    containing = expand_obstacle(&(in_area[0]), i, bbox, type);
    if (containing) {
      return (0);
    }
  }
  out_n = 0;
  for (i = 0; i < in_n; ++i) {
    tmp_area[0] = in_area[i];
    modify_area_sub0(tmp_area, &tmp_n);
    append_area_array(tmp_area, tmp_n, out_area, &out_n);
    assert(out_n <= N_TempArea);
  }
  return (out_n);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_rec_JS_line(AreaType* area)
{
  int linetype;
  double bbox[4];

  get_mr_bbox(bbox, area);
  linetype = calc_line_type(area->line[0][0], area->line[0][1]);
  if (linetype == VERTICAL) {
    assert(area->line[0][0].y > area->line[0][1].y);
    assert(area->line[1][0].y > area->line[1][1].y);
    area->line[0][0].y = area->line[1][0].y = bbox[NN];
    area->line[0][1].y = area->line[1][1].y = bbox[SS];
  } else if (linetype == HORIZONTAL) {
    assert(area->line[0][0].x > area->line[0][1].x);
    assert(area->line[1][0].x > area->line[1][1].x);
    area->line[0][0].x = area->line[1][0].x = bbox[EE];
    area->line[0][1].x = area->line[1][1].x = bbox[WW];
  } else {
    assert(0 == 1);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void check_rectilinear_line(AreaType* area)
{
  int linetype, linetype2;

  if (same_Point(area->line[0][0], area->line[0][1])) {
    assert(same_Point(area->line[1][0], area->line[1][1]));
    assert(equal(area->line[0][0].x, area->line[1][0].x)
           || equal(area->line[0][0].y, area->line[1][0].y));
  } else {
    linetype = calc_line_type(area->line[0][0], area->line[0][1]);
    linetype2 = calc_line_type(area->line[1][0], area->line[1][1]);
    assert(linetype == VERTICAL || linetype == HORIZONTAL);
    assert(linetype == linetype2);

    if (linetype == VERTICAL) {
      assert(area->line[0][0].y > area->line[0][1].y);
      assert(area->line[1][0].y > area->line[1][1].y);
    } else {
      assert(area->line[0][0].x > area->line[0][1].x);
      assert(area->line[1][0].x > area->line[1][1].x);
    }
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_rectilinear_line(AreaType* area, int linetype, int i, double z)
{
  if (linetype == VERTICAL) {
    area->line[0][i].y = area->line[1][i].y = z;
  } else {
    area->line[0][i].x = area->line[1][i].x = z;
  }
  check_rectilinear_line(area);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_area_case2_sub3(AreaType* area,
                            AreaType area_out[],
                            int* n_out,
                            double bbox[],
                            double obox[],
                            int d1,
                            int d2,
                            int linetype)
{
  area_out[*n_out] = *area;
  check_rectilinear_line(area);

  if (obox[d1] < bbox[d1] + FUZZ && obox[d2] > bbox[d2] - FUZZ) {
    modify_rectilinear_line(&(area_out[*n_out]), linetype, 1, obox[d1]);
    (*n_out)++;
    area_out[*n_out] = *area;
    modify_rectilinear_line(&(area_out[*n_out]), linetype, 0, obox[d2]);
    (*n_out)++;
  } else if (obox[d1] < bbox[d1] + FUZZ) {
    assert(obox[d2] <= bbox[d2] - FUZZ);
    modify_rectilinear_line(&(area_out[*n_out]), linetype, 1, obox[d1]);
    (*n_out)++;
  } else if (obox[d2] > bbox[d2] - FUZZ) {
    assert(obox[d1] >= bbox[d1] + FUZZ);
    modify_rectilinear_line(&(area_out[*n_out]), linetype, 0, obox[d2]);
    (*n_out)++;
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_area_case2_sub2(AreaType* area,
                            AreaType area_out[],
                            int* n_out,
                            double bbox[],
                            double obox[])
{
  int linetype, d1, d2;

  linetype = check_area_line(area);

  if (linetype == VERTICAL) {
    d1 = NN;
    d2 = SS;
  } else {
    d1 = EE;
    d2 = WW;
    assert(linetype == HORIZONTAL);
  }
  if (obox[d1] >= bbox[d1] + FUZZ && obox[d2] <= bbox[d2] - FUZZ) {
  } else {
    modify_area_case2_sub3(area, area_out, n_out, bbox, obox, d1, d2, linetype);
  }
  if (0) {
    check_area_array(area_out, *n_out);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int modify_area_case2_sub(AreaType* area,
                          AreaType area_out[],
                          int* n_out,
                          double obox[])
{
  int intersect;
  double bbox[4];

  get_JR_bbox(bbox, area);

  intersect = _bbox_overlap(bbox[WW],
                            bbox[SS],
                            bbox[EE],
                            bbox[NN],
                            obox[WW],
                            obox[SS],
                            obox[EE],
                            obox[NN],
                            0);
  if (intersect == NO) {
    check_area_line(area);
    area_out[*n_out] = *area;
    check_area_line(&(area_out[*n_out]));
    (*n_out)++;
  } else {
    modify_area_case2_sub2(area, area_out, n_out, bbox, obox);
  }
  return (intersect);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_mr_case2_sub(AreaType* area)
{
  int i, intersect;
  double bbox[4];

  get_JR_bbox(bbox, area);
  for (i = 0; i < N_Obstacle; ++i) {
    /*
        get_block_bbox(obox, &(Obstacle[i]));
    */
    intersect = _bbox_overlap(bbox[WW],
                              bbox[SS],
                              bbox[EE],
                              bbox[NN],
                              ObsBox[i][WW],
                              ObsBox[i][SS],
                              ObsBox[i][EE],
                              ObsBox[i][NN],
                              -FUZZ);
    assert(intersect == NO);
  }
  JS[0][0] = area->line[0][0];
  JS[0][1] = area->line[0][1];
  JS[1][0] = area->line[1][0];
  JS[1][1] = area->line[1][1];
  calc_JS_delay(area, area->area_L, area->area_R);

  if (Manhattan_arc(JS[0][0], JS[0][1])) {
    assert(same_Point(JS[0][0], JS[0][1]));
    assert(same_Point(JS[1][0], JS[1][1]));
    L_MS.MakeDiamond(JS[0][0], 0);
    R_MS.MakeDiamond(JS[1][0], 0);
  }
  MergeArea(area, area->area_L, area->area_R, NORMAL);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_mr_case2(int n, AreaType area[])
{
  int i;

  for (i = 0; i < n; ++i) {
    modify_mr_case2_sub(&(area[i]));
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int modify_area_case2(AreaType* area, AreaType area_out[])
{
  int n_in, n_out, i, j, change = NO, intersect;
  AreaType area_in[1000];

  area_out[0] = *area;
  n_out = 1;
  modify_rec_JS_line(&(area_out[0]));

  for (i = 0; i < N_Obstacle; ++i) {
    n_in = cp_area_array(area_out, n_out, area_in);
    if (0) {
      check_area_array(area_out, n_out);
    }
    if (0) {
      check_area_array(area_in, n_in);
    }
    n_out = 0;
    /*
        block = Obstacle[i];
        get_block_bbox(obox, &block);
    */
    for (j = 0; j < n_in; ++j) {
      intersect
          = modify_area_case2_sub(&(area_in[j]), area_out, &n_out, ObsBox[i]);
      assert(n_out < 50);
      if (intersect == YES) {
        change = YES;
      }
    }
    if (n_out == 0) {
      return (0);
    }
    assert(n_out < 50);
  }
  if (change) {
    modify_mr_case2(n_out, area_out);
  } else {
    assert(n_out == 1);
  }
  return (n_out);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int calc_path_balance_pt(AreaType* area,
                         PointType path[],
                         int n,
                         double delay[2][100],
                         double cap[2][100])
{
  double x, r, c;
  int i;

  delay[0][0] = path[0].max;
  delay[1][n - 1] = path[n - 1].max;
  cap[0][0] = area->area_L->capac;
  cap[1][n - 1] = area->area_R->capac;
  r = PURES[H_];
  c = PUCAP[H_];

  for (i = 1; i < n; ++i) {
    x = Point_dist(path[i], path[i - 1]);
    delay[0][i] = delay[0][i - 1] + r * x * (c * x / 2 + cap[0][i - 1]);
    cap[0][i] = cap[0][i - 1] + c * x;
  }
  for (i = n - 2; i >= 0; --i) {
    x = Point_dist(path[i], path[i + 1]);
    delay[1][i] = delay[1][i + 1] + r * x * (c * x / 2 + cap[1][i + 1]);
    cap[1][i] = cap[1][i + 1] + c * x;
  }
  for (i = 0; i < n - 2; ++i) {
    if (delay[0][i + 1] >= delay[1][i + 1]) {
      break;
    }
  }
  return (i);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void calc_new_path_delay(AreaType* area, PointType path[], int n)
{
  int i;
  PointType p0, p1;
  double x0, x1, t0, t1, r, c;

  r = PURES[H_];
  c = PUCAP[H_];

  p0 = path[0];
  p1 = path[n - 1];
  for (i = 0; i < n; ++i) {
    x0 = calc_subpath_length(path, 0, i, MANHATTAN);
    x1 = calc_subpath_length(path, i, n - 1, MANHATTAN);
    t0 = r * x0 * (c * x0 / 2 + area->area_L->capac);
    t1 = r * x1 * (c * x1 / 2 + area->area_R->capac);
    path[i].max = max(p0.max + t0, p1.max + t1);
    path[i].min = min(p0.min + t0, p1.min + t1);
  }
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void calc_feasible_subpath(PointType path[], int n, int k, int m[2])
{
  int i, j;
  double leng, dist;

  i = k;
  j = k + 1;
  leng = dist = Point_dist(path[k], path[k + 1]);
  while (i > 0) {
    dist = Point_dist(path[i - 1], path[j]);
    leng += Point_dist(path[i - 1], path[i]);
    if (!equal(leng, dist)) {
      break;
    }
    i--;
  }
  while (j < n - 1) {
    dist = Point_dist(path[i], path[j + 1]);
    leng += Point_dist(path[j], path[j + 1]);
    if (!equal(leng, dist)) {
      break;
    }
    j++;
  }
  m[0] = i;
  m[1] = j;
  assert(i < j);
  assert(i >= 0);
  assert(j < n);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void new_area_by_detour_sub(AreaType* area,
                            PointType path[],
                            int n,
                            int k,
                            double delay[2][100],
                            double cap[2][100])
{
  double r, c, d0, d1, t0, t1;
  int i;

  double x = Point_dist(path[k], path[k + 1]);
  calc_merge_distance(r,
                      c,
                      cap[0][k],
                      delay[0][k],
                      cap[1][k + 1],
                      delay[1][k + 1],
                      x,
                      &d0,
                      &d1);
  t0 = delay[0][k] + r * d0 * (c * d0 / 2 + cap[0][k]);
  t1 = delay[1][k + 1] + r * d1 * (c * d1 / 2 + cap[1][k + 1]);
  assert(equal(t0, t1));
  if (equal(d0, 0)) {
    area->mr[1] = area->mr[0] = path[k];
  } else if (equal(d1, 0)) {
    area->mr[1] = area->mr[0] = path[k + 1];
  } else {
    area->mr[0].x = (path[k].x * d1 + path[k + 1].x * d0) / (d0 + d1);
    area->mr[0].y = (path[k].y * d1 + path[k + 1].y * d0) / (d0 + d1);
    area->mr[1] = area->mr[0];
  }

  area->L_EdgeLen = d0;
  for (i = 0; i < k; ++i) {
    area->L_EdgeLen += Point_dist(path[i], path[i + 1]);
  }
  area->R_EdgeLen = d1;
  for (i = k + 1; i < n - 1; ++i) {
    area->R_EdgeLen += Point_dist(path[i], path[i + 1]);
  }

  x = area->L_EdgeLen;
  t0 = r * x * (c * x / 2 + area->area_L->capac);
  x = area->R_EdgeLen;
  t1 = r * x * (c * x / 2 + area->area_R->capac);
  assert(equal(path[0].max + t0, path[n - 1].max + t1));

  calc_n_vertex(area);
  for (i = 0; i < area->n_mr; ++i) {
    area->mr[i].max = path[0].max + t0;
    area->mr[i].min = min(path[0].min + t0, path[n - 1].min + t1);
    assert(pt_skew(area->mr[i]) <= Skew_B + FUZZ);
  }
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int new_area_by_detour_sub2(AreaType* area,
                            PointType path[],
                            int n,
                            int k,
                            AreaType out_area[])
{
  int i, m[2], n_area;
  double x0, x1, t0, t1;
  PointType pt0, pt1;
  AreaType tmp_area, tmp_area_L, tmp_area_R;

  calc_feasible_subpath(path, n, k, m);
  pt0 = path[m[0]];
  pt1 = path[m[1]];
  x0 = calc_subpath_length(path, 0, m[0], MANHATTAN);
  x1 = calc_subpath_length(path, m[1], n - 1, MANHATTAN);
  t0 = PURES[H_] * x0 * (PUCAP[H_] * x0 / 2 + area->area_L->capac);
  t1 = PURES[H_] * x1 * (PUCAP[H_] * x1 / 2 + area->area_R->capac);
  pt0.max = path[0].max + t0;
  pt0.min = path[0].min + t0;
  pt1.max = path[n - 1].max + t1;
  pt1.min = path[n - 1].min + t1;
  tmp_area_L = tmp_area_R = Node[0].area[0];
  tmp_area_L.capac = area->area_L->capac + x0 * PUCAP[H_];
  tmp_area_R.capac = area->area_R->capac + x1 * PUCAP[H_];
  tmp_area_L.mr[0] = pt0;
  tmp_area_R.mr[0] = pt1;
  JS[0][0] = JS[0][1] = pt0;
  JS[1][0] = JS[1][1] = pt1;
  L_MS.MakeDiamond(JS[0][0], 0);
  R_MS.MakeDiamond(JS[1][0], 0);

  MergeArea(&tmp_area, &tmp_area_L, &tmp_area_R, NORMAL);

  n_area = modify_blocked_areas_no_detour(&tmp_area, out_area);
  for (i = 0; i < n_area; ++i) {
    out_area[i].line[0][1] = out_area[i].line[0][0] = area->line[0][0];
    out_area[i].line[1][1] = out_area[i].line[1][0] = area->line[1][0];
    out_area[i].area_L = area->area_L;
    out_area[i].area_R = area->area_R;
    if (out_area[i].L_EdgeLen >= 0) {
      assert(out_area[i].R_EdgeLen >= 0);
      out_area[i].L_EdgeLen += x0;
      out_area[i].R_EdgeLen += x1;
    }
    MergeArea_sub(&(out_area[i]));
  }
  return (n_area);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int modify_blocked_areas_by_detour(AreaType* area, AreaType out_area[])
{
  int i, n, n_area = 1;
  PointType path[100], p1, p2, pt[2];
  double delay[2][100], cap[2][100];
  double cost;

  linedist(area->line[0][0],
           area->line[0][1],
           area->line[1][0],
           area->line[1][1],
           pt);
  cost = path_between_JSline_sub(area, area->line, path, &n, pt);
  assert(n > 2 && n <= 100);

  calc_BS_located(&(path[0]), area->area_L, &p1, &p2);
  calc_pt_delays(area->area_L, &(path[0]), p1, p2);

  calc_BS_located(&(path[n - 1]), area->area_R, &p1, &p2);
  calc_pt_delays(area->area_R, &(path[n - 1]), p1, p2);

  area->line[0][1] = area->line[0][0] = path[0];
  area->line[1][1] = area->line[1][0] = path[n - 1];

  i = calc_path_balance_pt(area, path, n, delay, cap);
  new_area_by_detour_sub(area, path, n, i, delay, cap);
  MergeArea_sub(area);
  out_area[0] = *area;
  /*
    n_area = new_area_by_detour_sub2(area, path, n,i, out_area);
  */
  return (n_area);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int modify_blocked_areas_no_detour(AreaType* area, AreaType new_area[])
{
  int linetype, n;

  linetype = calc_line_type(area->line[0][0], area->line[0][1]);
  if (linetype == MANHATTAN_ARC) {
    if (equal(area->dist, 0)) {
      n = 1;
      new_area[0] = *area;
    } else {
      n = modify_area_case1(area, new_area);
    }
  } else {
    assert(linetype == VERTICAL || linetype == HORIZONTAL);
    n = modify_area_case2(area, new_area);
  }

  return (n);
}

/***********************************************************************/
/*                                                                     */
/***********************************************************************/
int modify_blocked_areas_sub(AreaType area[], int n, AreaType area2[], int type)
{
  int i, n1, n2;
  AreaType area1[1000];

  n2 = 0;
  for (i = 0; i < n; ++i) {
    if (type == 0) {
      n1 = modify_blocked_areas_no_detour(&(area[i]), area1);
      if (0 && n1 == 0) {
        n1 = modify_blocked_areas_by_detour(&(area[i]), area1);
      }
    } else {
      n1 = modify_blocked_areas_by_detour(&(area[i]), area1);
    }
    assert(n1 < 1000);
    append_area_array(area1, n1, area2, &n2);
    assert(n2 < 1000);
  }
  return (n2);
}
/***********************************************************************/
/*                                                                     */
/***********************************************************************/
void modify_blocked_areas(AreaType area[], int* n, int b1, int b2)
{
  AreaType tmp_area[1000];

  if (N_Obstacle == 0) {
    return;
  }

  Buffered[0] = b1;
  Buffered[1] = b2;

  int n2 = modify_blocked_areas_sub(area, *n, tmp_area, 0);
  if (n2 == 0) {
    *n = min(*n, N_Area_Per_Node);
    n2 = modify_blocked_areas_sub(area, *n, tmp_area, 1);
  }
  *n = cp_area_array(tmp_area, n2, area);
}