bug #901655 ST_BUFFER asserts with a coplicated shape.

        Coinciding nodes can appear as a result of DOUBLE inaccuracy.
        We should test that before we start the loop.

        Also the spatial relations can be calculated faster if we check
        MBR relations first. And we do have the shape's MBR-s now.

per-file comments:
  sql/gcalc_slicescan.cc
        set_extent() method added.
bug #901655 ST_BUFFER asserts with a coplicated shape.
  sql/gcalc_slicescan.h
        set_extent() method declared.
bug #901655 ST_BUFFER asserts with a coplicated shape.
  sql/gcalc_tools.cc
bug #901655 ST_BUFFER asserts with a coplicated shape.
        checks for equal nodes added.
  sql/item_geofunc.cc
bug #901655 ST_BUFFER asserts with a coplicated shape.
        MBR for the shapes calculated, and MBR checks added before we
        start the heavy calculations.
  sql/spatial.h
bug #901655 ST_BUFFER asserts with a coplicated shape.
        MBR::buffer() method implemented.

5 files changed, 132 insertions, 192 deletions

diff --git a/sql/gcalc_slicescan.cc b/sql/gcalc_slicescan.cc
index 7021902ffb7..3c02bf62e53 100644
--- a/sql/gcalc_slicescan.cc
+++ b/sql/gcalc_slicescan.cc
@@ -406,7 +406,7 @@ static int de_weak_check(long double a, long double b, long double ex)
 
 static int de_check(long double a, long double b)
 {
-  return de_weak_check(a, b, (long double) 1e-10);
+  return de_weak_check(a, b, (long double) 1e-9);
 }
 #endif /*GCALC_CHECK_WITH_FLOAT*/
 
@@ -792,10 +792,51 @@ static int cmp_intersections(const Gcalc_heap::Info *i1,
 /* Internal coordinates implementation end */
 
 
+#define GCALC_SCALE_1 1e18
+
+static double find_scale(double extent)
+{
+  double scale= 1e-2;
+  while (scale < extent)
+    scale*= (double ) 10;
+  return GCALC_SCALE_1 / scale / 10;
+}
+
+
+void Gcalc_heap::set_extent(double xmin, double xmax, double ymin, double ymax)
+{
+  xmin= fabs(xmin);
+  xmax= fabs(xmax);
+  ymin= fabs(ymin);
+  ymax= fabs(ymax);
+
+  if (xmax < xmin)
+    xmax= xmin;
+  if (ymax < ymin)
+    ymax= ymin;
+
+  coord_extent= xmax > ymax ? xmax : ymax;
+  coord_extent= find_scale(coord_extent);
+#ifdef GCALC_CHECK_WITH_FLOAT
+  gcalc_coord_extent= &coord_extent;
+#endif /*GCALC_CHECK_WITH_FLOAT*/
+}
+
+
+void Gcalc_heap::free_point_info(Gcalc_heap::Info *i,
+                                 Gcalc_dyn_list::Item **i_hook)
+{
+  if (m_hook == &i->next)
+    m_hook= i_hook;
+  *i_hook= i->next;
+  free_item(i);
+  m_n_points--;
+}
+
+
 Gcalc_heap::Info *Gcalc_heap::new_point_info(double x, double y,
                                              gcalc_shape_info shape)
 {
-  double abs= fabs(x);
   Info *result= (Info *)new_item();
   if (!result)
     return NULL;
@@ -806,17 +847,8 @@ Gcalc_heap::Info *Gcalc_heap::new_point_info(double x, double y,
   result->shape= shape;
   result->top_node= 1;
   result->type= nt_shape_node;
-  if (m_n_points)
-  {
-    if (abs > coord_extent)
-      coord_extent= abs;
-  }
-  else
-    coord_extent= abs;
-
-  abs= fabs(y);
-  if (abs > coord_extent)
-    coord_extent= abs;
+  gcalc_set_double(result->ix, x, coord_extent);
+  gcalc_set_double(result->iy, y, coord_extent);
 
   m_n_points++;
   return result;
@@ -927,32 +959,12 @@ static int compare_point_info(const void *e0, const void *e1)
 }
 
 
-#define GCALC_SCALE_1 1e18
-
-static double find_scale(double extent)
-{
-  double scale= 1e-2;
-  while (scale < extent)
-    scale*= (double ) 10;
-  return GCALC_SCALE_1 / scale / 10;
-}
-
-
 void Gcalc_heap::prepare_operation()
 {
   Info *cur;
   GCALC_DBUG_ASSERT(m_hook);
-  coord_extent= find_scale(coord_extent);
-#ifdef GCALC_CHECK_WITH_FLOAT
-  gcalc_coord_extent= &coord_extent;
-#endif /*GCALC_CHECK_WITH_FLOAT*/
   *m_hook= NULL;
   m_hook= NULL; /* just to check it's not called twice */
-  for (cur= get_first(); cur; cur= cur->get_next())
-  {
-    gcalc_set_double(cur->ix, cur->x, coord_extent);
-    gcalc_set_double(cur->iy, cur->y, coord_extent);
-  }
   m_first= sort_list(compare_point_info, m_first, m_n_points);
 
   /* TODO - move this to the 'normal_scan' loop */
@@ -990,18 +1002,28 @@ int Gcalc_shape_transporter::int_add_point(gcalc_shape_info Info,
                                            double x, double y)
 {
   Gcalc_heap::Info *point;
-  GCALC_DBUG_ASSERT(!m_prev || m_prev->x != x || m_prev->y != y);
+  Gcalc_dyn_list::Item **hook;
+
+  hook= m_heap->get_cur_hook();
 
   if (!(point= m_heap->new_point_info(x, y, Info)))
     return 1;
   if (m_first)
   {
+    if (cmp_point_info(m_prev, point) == 0)
+    {
+      /* Coinciding points, do nothing */
+      m_heap->free_point_info(point, hook);
+      return 0;
+    }
+    GCALC_DBUG_ASSERT(!m_prev || m_prev->x != x || m_prev->y != y);
     m_prev->left= point;
     point->right= m_prev;
   }
   else
     m_first= point;
   m_prev= point;
+  m_prev_hook= hook;
   return 0;
 }
 
@@ -1029,10 +1051,20 @@ void Gcalc_shape_transporter::int_complete()
     return;
   }
 
-  GCALC_DBUG_ASSERT(m_prev->x != m_first->x || m_prev->y != m_first->y);
   /* polygon */
-  m_first->right= m_prev;
-  m_prev->left= m_first;
+  if (cmp_point_info(m_first, m_prev) == 0)
+  {
+    /* Coinciding points, remove the last one from the list */
+    m_prev->right->left= m_first;
+    m_first->right= m_prev->right;
+    m_heap->free_point_info(m_prev, m_prev_hook);
+  }
+  else
+  {
+    GCALC_DBUG_ASSERT(m_prev->x != m_first->x || m_prev->y != m_first->y);
+    m_first->right= m_prev;
+    m_prev->left= m_first;
+  }
 }
 
 
diff --git a/sql/gcalc_slicescan.h b/sql/gcalc_slicescan.h
index 9fc4cea5199..55de497f1ee 100644
--- a/sql/gcalc_slicescan.h
+++ b/sql/gcalc_slicescan.h
@@ -229,7 +229,9 @@ public:
     Gcalc_dyn_list(blk_size, sizeof(Info)),
     m_hook(&m_first), m_n_points(0)
   {}
+  void set_extent(double xmin, double xmax, double ymin, double ymax);
   Info *new_point_info(double x, double y, gcalc_shape_info shape);
+  void free_point_info(Info *i, Gcalc_dyn_list::Item **i_hook);
   Info *new_intersection(const Info *p1, const Info *p2,
                          const Info *p3, const Info *p4);
   void prepare_operation();
@@ -242,6 +244,8 @@ public:
   long double get_double(const Gcalc_internal_coord *c) const;
 #endif /*GCALC_CHECK_WITH_FLOAT*/
   double coord_extent;
+  Gcalc_dyn_list::Item **get_cur_hook() { return m_hook; }
+
 private:
   Gcalc_dyn_list::Item *m_first;
   Gcalc_dyn_list::Item **m_hook;
@@ -269,6 +273,7 @@ class Gcalc_shape_transporter
 private:
   Gcalc_heap::Info *m_first;
   Gcalc_heap::Info *m_prev;
+  Gcalc_dyn_list::Item **m_prev_hook;
   int m_shape_started;
   void int_complete();
 protected:
diff --git a/sql/gcalc_tools.cc b/sql/gcalc_tools.cc
index 161fa385f16..8af94039c2d 100644
--- a/sql/gcalc_tools.cc
+++ b/sql/gcalc_tools.cc
@@ -650,15 +650,6 @@ Gcalc_operation_reducer(Gcalc_function *fn, modes mode, size_t blk_size) :
 }
 
 
-#ifdef TMP_BLOCK
-void Gcalc_operation_reducer::res_point::set(const Gcalc_scan_iterator *si)
-{
-  if ((intersection_point= si->intersection_step()))
-    ii= si->get_cur_ii();
-  else
-    pi= si->get_cur_pi();
-}
-#endif /*TMP_BLOCK*/
 void Gcalc_operation_reducer::res_point::set(const Gcalc_scan_iterator *si)
 {
   intersection_point= si->intersection_step();
diff --git a/sql/item_geofunc.cc b/sql/item_geofunc.cc
index acaa50536ff..0c9ebc06749 100644
--- a/sql/item_geofunc.cc
+++ b/sql/item_geofunc.cc
@@ -657,140 +657,6 @@ static double distance_points(const Gcalc_heap::Info *a,
 }
 
 
-/*
-  Calculates the distance between objects.
-*/
-
-#ifdef TMP_BLOCK
-static int calc_distance(double *result, Gcalc_heap *collector, uint obj2_si,
-                         Gcalc_function *func, Gcalc_scan_iterator *scan_it)
-{
-  bool above_cur_point, cur_point_edge;
-  const Gcalc_scan_iterator::point *evpos;
-  const Gcalc_heap::Info *cur_point= NULL;
-  const Gcalc_heap::Info *dist_point;
-  const Gcalc_scan_iterator::point *ev;
-  double t, distance, cur_distance;
-  double ex, ey, vx, vy, e_sqrlen;
-  int o1, o2;
-
-  DBUG_ENTER("calc_distance");
-
-  above_cur_point= false;
-  distance= DBL_MAX;
-
-  while (scan_it->more_points())
-  {
-    if (scan_it->step())
-      goto mem_error;
-    evpos= scan_it->get_event_position();
-    ev= scan_it->get_events();
-    cur_point= NULL;
-
-    if (ev->simple_event())
-    {
-      cur_point= ev->pi;
-      goto calculate_distance;
-    }
-
-    /*
-       handling intersection we only need to check if it's the intersecion
-       of objects 1 and 2. In this case distance is 0
-    */
-    o1= 0;
-    o2= 0;
-    for (; ev; ev= ev->get_next())
-    {
-      if (ev->event != scev_intersection)
-        cur_point= ev->pi;
-      if (ev->pi->shape >= obj2_si)
-        o2= 1;
-      else
-        o1= 1;
-      if (o1 && o2)
-      {
-        distance= 0;
-        goto exit;
-      }
-    }
-    if (!cur_point)
-      continue;
-
-#ifdef TO_REMOVE
-    goto calculate_distance;
-    /*
-       having these events we need to check for possible intersection
-       of objects
-       scev_thread | scev_two_threads | scev_single_point
-    */
-    DBUG_ASSERT(ev & (scev_thread | scev_two_threads | scev_single_point));
-
-    func->clear_state();
-    for (Gcalc_point_iterator pit(scan_it); pit.point() != evpos; ++pit)
-    {
-      gcalc_shape_info si= pit.point()->get_shape();
-      if ((func->get_shape_kind(si) == Gcalc_function::shape_polygon))
-        func->invert_state(si);
-    }
-    func->invert_state(evpos->get_shape());
-    if (func->count())
-    {
-      /* Point of one object is inside the other - intersection found */
-      distance= 0;
-      goto exit;
-    }
-#endif /*TO_REMOVE*/
-
-calculate_distance:
-    if (cur_point->shape >= obj2_si)
-      continue;
-    cur_point_edge= !cur_point->is_bottom();
-
-    for (dist_point= collector->get_first(); dist_point; dist_point= dist_point->get_next())
-    {
-      /* We only check vertices of object 2 */
-      if (dist_point->shape < obj2_si)
-        continue;
-
-      /* if we have an edge to check */
-      if (dist_point->left)
-      {
-        t= count_edge_t(dist_point, dist_point->left, cur_point,
-                        ex, ey, vx, vy, e_sqrlen);
-        if ((t > 0.0) && (t < 1.0))
-        {
-          cur_distance= distance_to_line(ex, ey, vx, vy, e_sqrlen);
-          if (distance > cur_distance)
-            distance= cur_distance;
-        }
-      }
-      if (cur_point_edge)
-      {
-        t= count_edge_t(cur_point, cur_point->left, dist_point,
-                        ex, ey, vx, vy, e_sqrlen);
-        if ((t > 0.0) && (t < 1.0))
-        {
-          cur_distance= distance_to_line(ex, ey, vx, vy, e_sqrlen);
-          if (distance > cur_distance)
-            distance= cur_distance;
-        }
-      }
-      cur_distance= distance_points(cur_point, dist_point);
-      if (distance > cur_distance)
-        distance= cur_distance;
-    }
-  }
-
-exit:
-  *result= distance;
-  DBUG_RETURN(0);
-
-mem_error:
-  DBUG_RETURN(1);
-}
-#endif /*TMP_BLOCK*/
-
-
 #define GIS_ZERO 0.00000000001
 
 longlong Item_func_spatial_rel::val_int()
@@ -804,6 +670,8 @@ longlong Item_func_spatial_rel::val_int()
   int result= 0;
   int mask= 0;
   uint shape_a, shape_b;
+  MBR umbr, mbr1, mbr2;
+  const char *c_end;
 
   res1= args[0]->val_str(&tmp_value1);
   res2= args[1]->val_str(&tmp_value2);
@@ -818,19 +686,35 @@ longlong Item_func_spatial_rel::val_int()
 	!(g2= Geometry::construct(&buffer2, res2->ptr(), res2->length())))))
     goto exit;
 
+  g1->get_mbr(&mbr1, &c_end);
+  g2->get_mbr(&mbr2, &c_end);
+
+  umbr= mbr1;
+  umbr.add_mbr(&mbr2);
+  collector.set_extent(umbr.xmin, umbr.xmax, umbr.ymin, umbr.ymax);
+
+  mbr1.buffer(1e-5);
+
   switch (spatial_rel) {
     case SP_CONTAINS_FUNC:
+      if (!mbr1.contains(&mbr2))
+        goto exit;
       mask= 1;
       func.add_operation(Gcalc_function::op_difference, 2);
       /* Mind the g2 goes first. */
       null_value= g2->store_shapes(&trn) || g1->store_shapes(&trn);
       break;
     case SP_WITHIN_FUNC:
+      mbr2.buffer(2e-5);
+      if (!mbr1.within(&mbr2))
+        goto exit;
       mask= 1;
       func.add_operation(Gcalc_function::op_difference, 2);
       null_value= g1->store_shapes(&trn) || g2->store_shapes(&trn);
       break;
     case SP_EQUALS_FUNC:
+      if (!mbr1.contains(&mbr2))
+        goto exit;
       mask= 1;
       func.add_operation(Gcalc_function::op_symdifference, 2);
       null_value= g1->store_shapes(&trn) || g2->store_shapes(&trn);
@@ -841,6 +725,8 @@ longlong Item_func_spatial_rel::val_int()
       null_value= g1->store_shapes(&trn) || g2->store_shapes(&trn);
       break;
     case SP_INTERSECTS_FUNC:
+      if (!mbr1.intersects(&mbr2))
+        goto exit;
       func.add_operation(Gcalc_function::op_intersection, 2);
       null_value= g1->store_shapes(&trn) || g2->store_shapes(&trn);
       break;
@@ -898,15 +784,6 @@ longlong Item_func_spatial_rel::val_int()
   scan_it.init(&collector);
   scan_it.killed= (int *) &(current_thd->killed);
 
-#ifdef TMP_BLOCK
-  if (spatial_rel == SP_EQUALS_FUNC)
-  {
-    result= (g1->get_class_info()->m_type_id == g1->get_class_info()->m_type_id) &&
-            func_equals();
-    goto exit;
-  }
-#endif /*TMP_BLOCK*/
-
   if (func.alloc_states())
     goto exit;
 
@@ -935,6 +812,8 @@ String *Item_func_spatial_operation::val_str(String *str_value)
   Geometry *g1, *g2;
   uint32 srid= 0;
   Gcalc_operation_transporter trn(&func, &collector);
+  MBR mbr1, mbr2;
+  const char *c_end;
 
   if (func.reserve_op_buffer(1))
     DBUG_RETURN(0);
@@ -943,14 +822,23 @@ String *Item_func_spatial_operation::val_str(String *str_value)
   if ((null_value=
        (args[0]->null_value || args[1]->null_value ||
 	!(g1= Geometry::construct(&buffer1, res1->ptr(), res1->length())) ||
-	!(g2= Geometry::construct(&buffer2, res2->ptr(), res2->length())) ||
-	g1->store_shapes(&trn) || g2->store_shapes(&trn))))
+	!(g2= Geometry::construct(&buffer2, res2->ptr(), res2->length())))))
   {
     str_value= 0;
     goto exit;
   }
 
+  g1->get_mbr(&mbr1, &c_end);
+  g2->get_mbr(&mbr2, &c_end);
+  mbr1.add_mbr(&mbr2);
+  collector.set_extent(mbr1.xmin, mbr1.xmax, mbr1.ymin, mbr1.ymax);
   
+  if ((null_value= g1->store_shapes(&trn) || g2->store_shapes(&trn)))
+  {
+    str_value= 0;
+    goto exit;
+  }
+
   collector.prepare_operation();
   if (func.alloc_states())
     goto exit;
@@ -1374,12 +1262,18 @@ String *Item_func_buffer::val_str(String *str_value)
   uint32 srid= 0;
   String *str_result= NULL;
   Transporter trn(&func, &collector, dist);
+  MBR mbr;
+  const char *c_end;
 
   null_value= 1;
   if (args[0]->null_value || args[1]->null_value ||
-      !(g= Geometry::construct(&buffer, obj->ptr(), obj->length())))
+      !(g= Geometry::construct(&buffer, obj->ptr(), obj->length())) ||
+      g->get_mbr(&mbr, &c_end))
     goto mem_error;
 
+  if (dist > 0.0)
+    mbr.buffer(dist);
+  collector.set_extent(mbr.xmin, mbr.xmax, mbr.ymin, mbr.ymax);
   /*
     If the distance given is 0, the Buffer function is in fact NOOP,
     so it's natural just to return the argument1.
@@ -1446,6 +1340,8 @@ longlong Item_func_issimple::val_int()
   Geometry *g;
   int result= 1;
   const Gcalc_scan_iterator::event_point *ev;
+  MBR mbr;
+  const char *c_end;
 
   DBUG_ENTER("Item_func_issimple::val_int");
   DBUG_ASSERT(fixed == 1);
@@ -1454,6 +1350,8 @@ longlong Item_func_issimple::val_int()
       !(g= Geometry::construct(&buffer, swkb->ptr(), swkb->length())))
     DBUG_RETURN(0);
 
+  g->get_mbr(&mbr, &c_end);
+  collector.set_extent(mbr.xmin, mbr.xmax, mbr.ymin, mbr.ymax);
 
   if (g->get_class_info()->m_type_id == Geometry::wkb_point)
     DBUG_RETURN(1);
@@ -1682,6 +1580,8 @@ double Item_func_distance::val_real()
   String *res2= args[1]->val_str(&tmp_value2);
   Geometry_buffer buffer1, buffer2;
   Geometry *g1, *g2;
+  MBR mbr1, mbr2;
+  const char *c_end;
 
 
   if ((null_value= (args[0]->null_value || args[1]->null_value ||
@@ -1689,6 +1589,11 @@ double Item_func_distance::val_real()
           !(g2= Geometry::construct(&buffer2, res2->ptr(), res2->length())))))
     goto mem_error;
 
+  g1->get_mbr(&mbr1, &c_end);
+  g2->get_mbr(&mbr2, &c_end);
+  mbr1.add_mbr(&mbr2);
+  collector.set_extent(mbr1.xmin, mbr1.xmax, mbr1.ymin, mbr1.ymax);
+
   if ((g1->get_class_info()->m_type_id == Geometry::wkb_point) &&
       (g2->get_class_info()->m_type_id == Geometry::wkb_point))
   {
diff --git a/sql/spatial.h b/sql/spatial.h
index d55fc639acc..3ca5ed9ae20 100644
--- a/sql/spatial.h
+++ b/sql/spatial.h
@@ -98,6 +98,13 @@ struct MBR
     if (mbr->ymax > ymax)
       ymax= mbr->ymax;
   }
+  void buffer(double d)
+  {
+    xmin-= d;
+    ymin-= d;
+    xmax+= d;
+    ymax+= d;
+  }
 
   int equals(const MBR *mbr)
   {