mirror/mariadb
1
0
Fork 0
mirror of https://github.com/MariaDB/server.git synced 2025-01-16 12:02:42 +01:00
Code Issues Projects Releases Packages Wiki Activity Actions
mariadb/sql/gcalc_slicescan.cc
Marko Mäkelä be85d3e61b Merge 10.2 into 10.3
2019-05-14 17:18:46 +03:00

2011 lines
52 KiB
C++
Raw Blame History

/* Copyright (c) 2000, 2010 Oracle and/or its affiliates. All rights reserved.
Copyright (C) 2011 Monty Program Ab.
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; version 2 of the License.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1335 USA */
#include "mariadb.h"
#include <my_sys.h>
#include <m_string.h>
#ifdef HAVE_SPATIAL
#include "gcalc_slicescan.h"
#define PH_DATA_OFFSET 8
#define coord_to_float(d) ((double) d)
#define coord_eq(a, b) (a == b)
typedef int (*sc_compare_func)(const void*, const void*);
#define LS_LIST_ITEM Gcalc_dyn_list::Item
#define LS_COMPARE_FUNC_DECL sc_compare_func compare,
#define LS_COMPARE_FUNC_CALL(list_el1, list_el2) (*compare)(list_el1, list_el2)
#define LS_NEXT(A) (A)->next
#define LS_SET_NEXT(A,val) (A)->next= val
#define LS_P_NEXT(A) &(A)->next
#define LS_NAME sort_list
#define LS_SCOPE static
#define LS_STRUCT_NAME sort_list_stack_struct
#include "plistsort.c"
#define GCALC_COORD_MINUS 0x80000000
#define FIRST_DIGIT(d) ((d) & 0x7FFFFFFF)
#define GCALC_SIGN(d) ((d) & 0x80000000)
static Gcalc_scan_iterator::point *eq_sp(const Gcalc_heap::Info *pi)
{
GCALC_DBUG_ASSERT(pi->type == Gcalc_heap::nt_eq_node);
return (Gcalc_scan_iterator::point *) pi->node.eq.data;
}
static Gcalc_scan_iterator::intersection_info *i_data(const Gcalc_heap::Info *pi)
{
GCALC_DBUG_ASSERT(pi->type == Gcalc_heap::nt_intersection);
return (Gcalc_scan_iterator::intersection_info *) pi->node.intersection.data;
}
#ifndef GCALC_DBUG_OFF
int gcalc_step_counter= 0;
void GCALC_DBUG_CHECK_COUNTER()
{
if (++gcalc_step_counter == 0)
GCALC_DBUG_PRINT(("step_counter_0"));
else
GCALC_DBUG_PRINT(("%d step_counter", gcalc_step_counter));
}
const char *gcalc_ev_name(int ev)
{
switch (ev)
{
case scev_none:
return "n";
case scev_thread:
return "t";
case scev_two_threads:
return "tt";
case scev_end:
return "e";
case scev_two_ends:
return "ee";
case scev_intersection:
return "i";
case scev_point:
return "p";
case scev_single_point:
return "sp";
default:;
};
GCALC_DBUG_ASSERT(0);
return "unk";
}
static int gcalc_pi_str(char *str, const Gcalc_heap::Info *pi, const char *postfix)
{
return sprintf(str, "%s %d %d | %s %d %d%s",
GCALC_SIGN(pi->node.shape.ix[0]) ? "-":"", FIRST_DIGIT(pi->node.shape.ix[0]),pi->node.shape.ix[1],
GCALC_SIGN(pi->node.shape.iy[0]) ? "-":"", FIRST_DIGIT(pi->node.shape.iy[0]),pi->node.shape.iy[1],
postfix);
}
static void GCALC_DBUG_PRINT_PI(const Gcalc_heap::Info *pi)
{
char buf[128];
int n_buf;
if (pi->type == Gcalc_heap::nt_intersection)
{
const Gcalc_scan_iterator::intersection_info *ic= i_data(pi);
GCALC_DBUG_PRINT(("intersection point %d %d",
ic->edge_a->thread, ic->edge_b->thread));
return;
}
if (pi->type == Gcalc_heap::nt_eq_node)
{
const Gcalc_scan_iterator::point *e= eq_sp(pi);
GCALC_DBUG_PRINT(("eq point %d", e->thread));
return;
}
n_buf= gcalc_pi_str(buf, pi, "");
buf[n_buf]= 0;
GCALC_DBUG_PRINT(("%s", buf));
}
static void GCALC_DBUG_PRINT_SLICE(const char *header,
const Gcalc_scan_iterator::point *slice)
{
size_t nbuf;
char buf[1024];
nbuf= strlen(header);
strcpy(buf, header);
for (; slice; slice= slice->get_next())
{
size_t lnbuf= nbuf;
lnbuf+= sprintf(buf + lnbuf, "%d\t", slice->thread);
lnbuf+= sprintf(buf + lnbuf, "%s\t", gcalc_ev_name(slice->event));
lnbuf+= gcalc_pi_str(buf + lnbuf, slice->pi, "\t");
if (slice->is_bottom())
lnbuf+= sprintf(buf+lnbuf, "bt\t");
else
lnbuf+= gcalc_pi_str(buf+lnbuf, slice->next_pi, "\t");
buf[lnbuf]= 0;
GCALC_DBUG_PRINT(("%s", buf));
}
}
#else
#define GCALC_DBUG_CHECK_COUNTER() do { } while(0)
#define GCALC_DBUG_PRINT_PI(pi) do { } while(0)
#define GCALC_DBUG_PRINT_SLICE(a, b) do { } while(0)
#define GCALC_DBUG_PRINT_INTERSECTIONS(a) do { } while(0)
#define GCALC_DBUG_PRINT_STATE(a) do { } while(0)
#endif /*GCALC_DBUG_OFF*/
Gcalc_dyn_list::Gcalc_dyn_list(size_t blk_size, size_t sizeof_item):
m_blk_size(blk_size - ALLOC_ROOT_MIN_BLOCK_SIZE),
m_sizeof_item(ALIGN_SIZE(sizeof_item)),
m_points_per_blk((uint)((m_blk_size - PH_DATA_OFFSET) / m_sizeof_item)),
m_blk_hook(&m_first_blk),
m_free(NULL),
m_keep(NULL)
{}
Gcalc_dyn_list::Gcalc_dyn_list(const Gcalc_dyn_list &dl)
{
m_blk_size= dl.m_blk_size;
m_sizeof_item= dl.m_sizeof_item;
m_points_per_blk= dl.m_points_per_blk;
m_blk_hook= &m_first_blk;
m_free= NULL;
m_keep= NULL;
}
void Gcalc_dyn_list::format_blk(void* block)
{
Item *pi_end, *cur_pi, *first_pi;
GCALC_DBUG_ASSERT(m_free == NULL);
first_pi= cur_pi= (Item *)(((char *)block) + PH_DATA_OFFSET);
pi_end= ptr_add(first_pi, m_points_per_blk - 1);
do {
cur_pi= cur_pi->next= ptr_add(cur_pi, 1);
} while (cur_pi<pi_end);
cur_pi->next= m_free;
m_free= first_pi;
}
Gcalc_dyn_list::Item *Gcalc_dyn_list::alloc_new_blk()
{
void *new_block= my_malloc(m_blk_size, MYF(MY_WME));
if (!new_block)
return NULL;
*m_blk_hook= new_block;
m_blk_hook= (void**)new_block;
format_blk(new_block);
return new_item();
}
static void free_blk_list(void *list)
{
void *next_blk;
while (list)
{
next_blk= *((void **)list);
my_free(list);
list= next_blk;
}
}
void Gcalc_dyn_list::cleanup()
{
*m_blk_hook= NULL;
free_blk_list(m_first_blk);
m_first_blk= NULL;
m_blk_hook= &m_first_blk;
m_free= NULL;
}
Gcalc_dyn_list::~Gcalc_dyn_list()
{
cleanup();
}
void Gcalc_dyn_list::reset()
{
*m_blk_hook= NULL;
if (m_first_blk)
{
free_blk_list(*((void **)m_first_blk));
m_blk_hook= (void**)m_first_blk;
m_free= NULL;
format_blk(m_first_blk);
}
}
/* Internal coordinate operations implementations */
void gcalc_set_zero(Gcalc_internal_coord *d, int d_len)
{
do
{
d[--d_len]= 0;
} while (d_len);
}
int gcalc_is_zero(const Gcalc_internal_coord *d, int d_len)
{
do
{
if (d[--d_len] != 0)
return 0;
} while (d_len);
return 1;
}
#ifdef GCALC_CHECK_WITH_FLOAT
static double *gcalc_coord_extent= NULL;
long double gcalc_get_double(const Gcalc_internal_coord *d, int d_len)
{
int n= 1;
long double res= (long double) FIRST_DIGIT(d[0]);
do
{
res*= (long double) GCALC_DIG_BASE;
res+= (long double) d[n];
} while(++n < d_len);
n= 0;
do
{
if ((n & 1) && gcalc_coord_extent)
res/= *gcalc_coord_extent;
} while(++n < d_len);
if (GCALC_SIGN(d[0]))
res*= -1.0;
return res;
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
static void do_add(Gcalc_internal_coord *result, int result_len,
const Gcalc_internal_coord *a,
const Gcalc_internal_coord *b)
{
int n_digit= result_len-1;
gcalc_digit_t carry= 0;
do
{
if ((result[n_digit]=
a[n_digit] + b[n_digit] + carry) >= GCALC_DIG_BASE)
{
carry= 1;
result[n_digit]-= GCALC_DIG_BASE;
}
else
carry= 0;
} while (--n_digit);
result[0]= (a[0] + FIRST_DIGIT(b[0]) + carry);
GCALC_DBUG_ASSERT(FIRST_DIGIT(result[0]) < GCALC_DIG_BASE);
}
static void do_sub(Gcalc_internal_coord *result, int result_len,
const Gcalc_internal_coord *a,
const Gcalc_internal_coord *b)
{
int n_digit= result_len-1;
gcalc_digit_t carry= 0;
gcalc_digit_t cur_b, cur_a;
do
{
cur_b= b[n_digit] + carry;
cur_a= a[n_digit];
if (cur_a < cur_b)
{
carry= 1;
result[n_digit]= (GCALC_DIG_BASE - cur_b) + cur_a;
}
else
{
carry= 0;
result[n_digit]= cur_a - cur_b;
}
} while (--n_digit);
result[0]= a[0] - FIRST_DIGIT(b[0]) - carry;
GCALC_DBUG_ASSERT(FIRST_DIGIT(a[0]) >= FIRST_DIGIT(b[0]) + carry);
GCALC_DBUG_ASSERT(!gcalc_is_zero(result, result_len));
}
/*
static void do_sub(Gcalc_internal_coord *result, int result_len,
const Gcalc_internal_coord *a,
const Gcalc_internal_coord *b)
{
int n_digit= result_len-1;
gcalc_digit_t carry= 0;
do
{
if ((result[n_digit]= a[n_digit] - b[n_digit] - carry) < 0)
{
carry= 1;
result[n_digit]+= GCALC_DIG_BASE;
}
else
carry= 0;
} while (--n_digit);
result[0]= a[0] - FIRST_DIGIT(b[0]) - carry;
GCALC_DBUG_ASSERT(FIRST_DIGIT(a[0]) - FIRST_DIGIT(b[0]) - carry >= 0);
GCALC_DBUG_ASSERT(!gcalc_is_zero(result, result_len));
}
*/
static int do_cmp(const Gcalc_internal_coord *a,
const Gcalc_internal_coord *b, int len)
{
int n_digit= 1;
if ((FIRST_DIGIT(a[0]) != FIRST_DIGIT(b[0])))
return FIRST_DIGIT(a[0]) > FIRST_DIGIT(b[0]) ? 1 : -1;
do
{
if ((a[n_digit] != b[n_digit]))
return a[n_digit] > b[n_digit] ? 1 : -1;
} while (++n_digit < len);
return 0;
}
#ifdef GCALC_CHECK_WITH_FLOAT
static int de_weak_check(long double a, long double b, long double ex)
{
long double d= a - b;
if (d < ex && d > -ex)
return 1;
d/= fabsl(a) + fabsl(b);
if (d < ex && d > -ex)
return 1;
return 0;
}
static int de_check(long double a, long double b)
{
return de_weak_check(a, b, (long double) 1e-9);
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
void gcalc_mul_coord(Gcalc_internal_coord *result, int result_len,
const Gcalc_internal_coord *a, int a_len,
const Gcalc_internal_coord *b, int b_len)
{
GCALC_DBUG_ASSERT(result_len == a_len + b_len);
GCALC_DBUG_ASSERT(a_len >= b_len);
int n_a, n_b, n_res;
gcalc_digit_t carry= 0;
gcalc_set_zero(result, result_len);
n_a= a_len - 1;
do
{
gcalc_coord2 cur_a= n_a ? a[n_a] : FIRST_DIGIT(a[0]);
n_b= b_len - 1;
do
{
gcalc_coord2 cur_b= n_b ? b[n_b] : FIRST_DIGIT(b[0]);
gcalc_coord2 mul= cur_a * cur_b + carry + result[n_a + n_b + 1];
result[n_a + n_b + 1]= mul % GCALC_DIG_BASE;
carry= (gcalc_digit_t) (mul / (gcalc_coord2) GCALC_DIG_BASE);
} while (n_b--);
if (carry)
{
for (n_res= n_a; (result[n_res]+= carry) >= GCALC_DIG_BASE;
n_res--)
{
result[n_res]-= GCALC_DIG_BASE;
carry= 1;
}
carry= 0;
}
} while (n_a--);
if (!gcalc_is_zero(result, result_len))
result[0]|= GCALC_SIGN(a[0] ^ b[0]);
#ifdef GCALC_CHECK_WITH_FLOAT
GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, a_len) *
gcalc_get_double(b, b_len),
gcalc_get_double(result, result_len)));
#endif /*GCALC_CHECK_WITH_FLOAT*/
}
inline void gcalc_mul_coord1(Gcalc_coord1 result,
const Gcalc_coord1 a, const Gcalc_coord1 b)
{
return gcalc_mul_coord(result, GCALC_COORD_BASE2,
a, GCALC_COORD_BASE, b, GCALC_COORD_BASE);
}
void gcalc_add_coord(Gcalc_internal_coord *result, int result_len,
const Gcalc_internal_coord *a,
const Gcalc_internal_coord *b)
{
if (GCALC_SIGN(a[0]) == GCALC_SIGN(b[0]))
do_add(result, result_len, a, b);
else
{
int cmp_res= do_cmp(a, b, result_len);
if (cmp_res == 0)
gcalc_set_zero(result, result_len);
else if (cmp_res > 0)
do_sub(result, result_len, a, b);
else
do_sub(result, result_len, b, a);
}
#ifdef GCALC_CHECK_WITH_FLOAT
GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, result_len) +
gcalc_get_double(b, result_len),
gcalc_get_double(result, result_len)));
#endif /*GCALC_CHECK_WITH_FLOAT*/
}
void gcalc_sub_coord(Gcalc_internal_coord *result, int result_len,
const Gcalc_internal_coord *a,
const Gcalc_internal_coord *b)
{
if (GCALC_SIGN(a[0] ^ b[0]))
do_add(result, result_len, a, b);
else
{
int cmp_res= do_cmp(a, b, result_len);
if (cmp_res == 0)
gcalc_set_zero(result, result_len);
else if (cmp_res > 0)
do_sub(result, result_len, a, b);
else
{
do_sub(result, result_len, b, a);
result[0]^= GCALC_COORD_MINUS;
}
}
#ifdef GCALC_CHECK_WITH_FLOAT
GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, result_len) -
gcalc_get_double(b, result_len),
gcalc_get_double(result, result_len)));
#endif /*GCALC_CHECK_WITH_FLOAT*/
}
inline void gcalc_sub_coord1(Gcalc_coord1 result,
const Gcalc_coord1 a, const Gcalc_coord1 b)
{
return gcalc_sub_coord(result, GCALC_COORD_BASE, a, b);
}
int gcalc_cmp_coord(const Gcalc_internal_coord *a,
const Gcalc_internal_coord *b, int len)
{
int n_digit= 0;
int result= 0;
do
{
if (a[n_digit] == b[n_digit])
{
n_digit++;
continue;
}
if (a[n_digit] > b[n_digit])
result= GCALC_SIGN(a[0]) ? -1 : 1;
else
result= GCALC_SIGN(b[0]) ? 1 : -1;
break;
} while (n_digit < len);
#ifdef GCALC_CHECK_WITH_FLOAT
if (result == 0)
GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len),
gcalc_get_double(b, len)));
else if (result == 1)
GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len),
gcalc_get_double(b, len)) ||
gcalc_get_double(a, len) > gcalc_get_double(b, len));
else
GCALC_DBUG_ASSERT(de_check(gcalc_get_double(a, len),
gcalc_get_double(b, len)) ||
gcalc_get_double(a, len) < gcalc_get_double(b, len));
#endif /*GCALC_CHECK_WITH_FLOAT*/
return result;
}
#define gcalc_cmp_coord1(a, b) gcalc_cmp_coord(a, b, GCALC_COORD_BASE)
int gcalc_set_double(Gcalc_internal_coord *c, double d, double ext)
{
int sign;
double ds= d * ext;
if ((sign= ds < 0))
ds= -ds;
c[0]= (gcalc_digit_t) (ds / (double) GCALC_DIG_BASE);
c[1]= (gcalc_digit_t) (ds - ((double) c[0]) * (double) GCALC_DIG_BASE);
if (c[1] >= GCALC_DIG_BASE)
{
c[1]= 0;
c[0]++;
}
if (sign && (c[0] | c[1]))
c[0]|= GCALC_COORD_MINUS;
#ifdef GCALC_CHECK_WITH_FLOAT
GCALC_DBUG_ASSERT(de_check(d, gcalc_get_double(c, 2)));
#endif /*GCALC_CHECK_WITH_FLOAT*/
return 0;
}
typedef gcalc_digit_t Gcalc_coord4[GCALC_COORD_BASE*4];
typedef gcalc_digit_t Gcalc_coord5[GCALC_COORD_BASE*5];
void Gcalc_scan_iterator::intersection_info::do_calc_t()
{
Gcalc_coord1 a2_a1x, a2_a1y;
Gcalc_coord2 x1y2, x2y1;
gcalc_sub_coord1(a2_a1x, edge_b->pi->node.shape.ix, edge_a->pi->node.shape.ix);
gcalc_sub_coord1(a2_a1y, edge_b->pi->node.shape.iy, edge_a->pi->node.shape.iy);
GCALC_DBUG_ASSERT(!gcalc_is_zero(edge_a->dy, GCALC_COORD_BASE) ||
!gcalc_is_zero(edge_b->dy, GCALC_COORD_BASE));
gcalc_mul_coord1(x1y2, edge_a->dx, edge_b->dy);
gcalc_mul_coord1(x2y1, edge_a->dy, edge_b->dx);
gcalc_sub_coord(t_b, GCALC_COORD_BASE2, x1y2, x2y1);
gcalc_mul_coord1(x1y2, a2_a1x, edge_b->dy);
gcalc_mul_coord1(x2y1, a2_a1y, edge_b->dx);
gcalc_sub_coord(t_a, GCALC_COORD_BASE2, x1y2, x2y1);
t_calculated= 1;
}
void Gcalc_scan_iterator::intersection_info::do_calc_y()
{
GCALC_DBUG_ASSERT(t_calculated);
Gcalc_coord3 a_tb, b_ta;
gcalc_mul_coord(a_tb, GCALC_COORD_BASE3,
t_b, GCALC_COORD_BASE2, edge_a->pi->node.shape.iy, GCALC_COORD_BASE);
gcalc_mul_coord(b_ta, GCALC_COORD_BASE3,
t_a, GCALC_COORD_BASE2, edge_a->dy, GCALC_COORD_BASE);
gcalc_add_coord(y_exp, GCALC_COORD_BASE3, a_tb, b_ta);
y_calculated= 1;
}
void Gcalc_scan_iterator::intersection_info::do_calc_x()
{
GCALC_DBUG_ASSERT(t_calculated);
Gcalc_coord3 a_tb, b_ta;
gcalc_mul_coord(a_tb, GCALC_COORD_BASE3,
t_b, GCALC_COORD_BASE2, edge_a->pi->node.shape.ix, GCALC_COORD_BASE);
gcalc_mul_coord(b_ta, GCALC_COORD_BASE3,
t_a, GCALC_COORD_BASE2, edge_a->dx, GCALC_COORD_BASE);
gcalc_add_coord(x_exp, GCALC_COORD_BASE3, a_tb, b_ta);
x_calculated= 1;
}
static int cmp_node_isc(const Gcalc_heap::Info *node,
const Gcalc_heap::Info *isc)
{
GCALC_DBUG_ASSERT(node->type == Gcalc_heap::nt_shape_node);
Gcalc_scan_iterator::intersection_info *inf= i_data(isc);
Gcalc_coord3 exp;
int result;
inf->calc_t();
inf->calc_y_exp();
gcalc_mul_coord(exp, GCALC_COORD_BASE3,
inf->t_b, GCALC_COORD_BASE2, node->node.shape.iy, GCALC_COORD_BASE);
result= gcalc_cmp_coord(exp, inf->y_exp, GCALC_COORD_BASE3);
#ifdef GCALC_CHECK_WITH_FLOAT
long double int_x, int_y;
isc->calc_xy_ld(&int_x, &int_y);
if (result < 0)
{
if (!de_check(int_y, node->node.shape.y) && node->node.shape.y > int_y)
GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g < %LG", node->node.shape.y, int_y));
}
else if (result > 0)
{
if (!de_check(int_y, node->node.shape.y) && node->node.shape.y < int_y)
GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g > %LG", node->node.shape.y, int_y));
}
else
{
if (!de_check(int_y, node->node.shape.y))
GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscy %g == %LG", node->node.shape.y, int_y));
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
if (result)
goto exit;
inf->calc_x_exp();
gcalc_mul_coord(exp, GCALC_COORD_BASE3,
inf->t_b, GCALC_COORD_BASE2, node->node.shape.ix, GCALC_COORD_BASE);
result= gcalc_cmp_coord(exp, inf->x_exp, GCALC_COORD_BASE3);
#ifdef GCALC_CHECK_WITH_FLOAT
if (result < 0)
{
if (!de_check(int_x, node->node.shape.x) && node->node.shape.x > int_x)
GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g < %LG",
node->node.shape.x, int_x));
}
else if (result > 0)
{
if (!de_check(int_x, node->node.shape.x) && node->node.shape.x < int_x)
GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g > %LG",
node->node.shape.x, int_x));
}
else
{
if (!de_check(int_x, node->node.shape.x))
GCALC_DBUG_PRINT(("floatcheck cmp_nod_iscx failed %g == %LG",
node->node.shape.x, int_x));
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
exit:
return result;
}
static int cmp_intersections(const Gcalc_heap::Info *i1,
const Gcalc_heap::Info *i2)
{
Gcalc_scan_iterator::intersection_info *inf1= i_data(i1);
Gcalc_scan_iterator::intersection_info *inf2= i_data(i2);
Gcalc_coord5 exp_a, exp_b;
int result;
inf1->calc_t();
inf2->calc_t();
inf1->calc_y_exp();
inf2->calc_y_exp();
gcalc_mul_coord(exp_a, GCALC_COORD_BASE5,
inf1->y_exp, GCALC_COORD_BASE3, inf2->t_b, GCALC_COORD_BASE2);
gcalc_mul_coord(exp_b, GCALC_COORD_BASE5,
inf2->y_exp, GCALC_COORD_BASE3, inf1->t_b, GCALC_COORD_BASE2);
result= gcalc_cmp_coord(exp_a, exp_b, GCALC_COORD_BASE5);
#ifdef GCALC_CHECK_WITH_FLOAT
long double x1, y1, x2, y2;
i1->calc_xy_ld(&x1, &y1);
i2->calc_xy_ld(&x2, &y2);
if (result < 0)
{
if (!de_check(y1, y2) && y2 > y1)
GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG < %LG",
y2, y1));
}
else if (result > 0)
{
if (!de_check(y1, y2) && y2 < y1)
GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG > %LG",
y2, y1));
}
else
{
if (!de_check(y1, y2))
GCALC_DBUG_PRINT(("floatcheck cmp_intersections_y failed %LG == %LG",
y2, y1));
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
if (result != 0)
return result;
inf1->calc_x_exp();
inf2->calc_x_exp();
gcalc_mul_coord(exp_a, GCALC_COORD_BASE5,
inf1->x_exp, GCALC_COORD_BASE3, inf2->t_b, GCALC_COORD_BASE2);
gcalc_mul_coord(exp_b, GCALC_COORD_BASE5,
inf2->x_exp, GCALC_COORD_BASE3, inf1->t_b, GCALC_COORD_BASE2);
result= gcalc_cmp_coord(exp_a, exp_b, GCALC_COORD_BASE5);
#ifdef GCALC_CHECK_WITH_FLOAT
if (result < 0)
{
if (!de_check(x1, x2) && x2 > x1)
GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG < %LG",
x2, x1));
}
else if (result > 0)
{
if (!de_check(x1, x2) && x2 < x1)
GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG > %LG",
x2, x1));
}
else
{
if (!de_check(x1, x2))
GCALC_DBUG_PRINT(("floatcheck cmp_intersectionsx failed %LG == %LG",
x2, x1));
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
return result;
}
/* 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)
{
Info *result= (Info *)new_item();
if (!result)
return NULL;
*m_hook= result;
m_hook= &result->next;
result->node.shape.x= x;
result->node.shape.y= y;
result->node.shape.shape= shape;
result->node.shape.top_node= 1;
result->type= nt_shape_node;
gcalc_set_double(result->node.shape.ix, x, coord_extent);
gcalc_set_double(result->node.shape.iy, y, coord_extent);
m_n_points++;
return result;
}
static Gcalc_heap::Info *new_intersection(
Gcalc_heap *heap, Gcalc_scan_iterator::intersection_info *ii)
{
Gcalc_heap::Info *isc= (Gcalc_heap::Info *)heap->new_item();
if (!isc)
return 0;
isc->type= Gcalc_heap::nt_intersection;
isc->node.intersection.p1= ii->edge_a->pi;
isc->node.intersection.p2= ii->edge_a->next_pi;
isc->node.intersection.p3= ii->edge_b->pi;
isc->node.intersection.p4= ii->edge_b->next_pi;
isc->node.intersection.data= ii;
return isc;
}
static Gcalc_heap::Info *new_eq_point(
Gcalc_heap *heap, const Gcalc_heap::Info *p,
Gcalc_scan_iterator::point *edge)
{
Gcalc_heap::Info *eqp= (Gcalc_heap::Info *)heap->new_item();
if (!eqp)
return 0;
eqp->type= Gcalc_heap::nt_eq_node;
eqp->node.eq.node= p;
eqp->node.eq.data= edge;
return eqp;
}
void Gcalc_heap::Info::calc_xy(double *x, double *y) const
{
double b0_x= node.intersection.p2->node.shape.x - node.intersection.p1->node.shape.x;
double b0_y= node.intersection.p2->node.shape.y - node.intersection.p1->node.shape.y;
double b1_x= node.intersection.p4->node.shape.x - node.intersection.p3->node.shape.x;
double b1_y= node.intersection.p4->node.shape.y - node.intersection.p3->node.shape.y;
double b0xb1= b0_x * b1_y - b0_y * b1_x;
double t= (node.intersection.p3->node.shape.x - node.intersection.p1->node.shape.x) * b1_y - (node.intersection.p3->node.shape.y - node.intersection.p1->node.shape.y) * b1_x;
t/= b0xb1;
*x= node.intersection.p1->node.shape.x + b0_x * t;
*y= node.intersection.p1->node.shape.y + b0_y * t;
}
#ifdef GCALC_CHECK_WITH_FLOAT
void Gcalc_heap::Info::calc_xy_ld(long double *x, long double *y) const
{
long double b0_x= ((long double) p2->node.shape.x) - p1->node.shape.x;
long double b0_y= ((long double) p2->node.shape.y) - p1->node.shape.y;
long double b1_x= ((long double) p4->node.shape.x) - p3->node.shape.x;
long double b1_y= ((long double) p4->node.shape.y) - p3->node.shape.y;
long double b0xb1= b0_x * b1_y - b0_y * b1_x;
long double ax= ((long double) p3->node.shape.x) - p1->node.shape.x;
long double ay= ((long double) p3->node.shape.y) - p1->node.shape.y;
long double t_a= ax * b1_y - ay * b1_x;
long double hx= (b0xb1 * (long double) p1->node.shape.x + b0_x * t_a);
long double hy= (b0xb1 * (long double) p1->node.shape.y + b0_y * t_a);
if (fabs(b0xb1) < 1e-15)
{
*x= p1->node.shape.x;
*y= p1->node.shape.y;
return;
}
*x= hx/b0xb1;
*y= hy/b0xb1;
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
static int cmp_point_info(const Gcalc_heap::Info *i0,
const Gcalc_heap::Info *i1)
{
int cmp_y= gcalc_cmp_coord1(i0->node.shape.iy, i1->node.shape.iy);
if (cmp_y)
return cmp_y;
return gcalc_cmp_coord1(i0->node.shape.ix, i1->node.shape.ix);
}
static inline void trim_node(Gcalc_heap::Info *node, Gcalc_heap::Info *prev_node)
{
if (!node)
return;
node->node.shape.top_node= 0;
GCALC_DBUG_ASSERT((node->node.shape.left == prev_node) || (node->node.shape.right == prev_node));
if (node->node.shape.left == prev_node)
node->node.shape.left= node->node.shape.right;
node->node.shape.right= NULL;
GCALC_DBUG_ASSERT(cmp_point_info(node, prev_node));
}
static int compare_point_info(const void *e0, const void *e1)
{
const Gcalc_heap::Info *i0= (const Gcalc_heap::Info *)e0;
const Gcalc_heap::Info *i1= (const Gcalc_heap::Info *)e1;
return cmp_point_info(i0, i1) > 0;
}
void Gcalc_heap::prepare_operation()
{
Info *cur;
GCALC_DBUG_ASSERT(m_hook);
*m_hook= NULL;
m_hook= NULL; /* just to check it's not called twice */
m_first= sort_list(compare_point_info, m_first, m_n_points);
/* TODO - move this to the 'normal_scan' loop */
for (cur= get_first(); cur; cur= cur->get_next())
{
trim_node(cur->node.shape.left, cur);
trim_node(cur->node.shape.right, cur);
}
}
void Gcalc_heap::reset()
{
if (m_n_points)
{
if (m_hook)
*m_hook= NULL;
free_list(m_first);
m_n_points= 0;
}
m_hook= &m_first;
}
int Gcalc_shape_transporter::int_single_point(gcalc_shape_info Info,
double x, double y)
{
Gcalc_heap::Info *point= m_heap->new_point_info(x, y, Info);
if (!point)
return 1;
point->node.shape.left= point->node.shape.right= 0;
return 0;
}
int Gcalc_shape_transporter::int_add_point(gcalc_shape_info Info,
double x, double y)
{
Gcalc_heap::Info *point;
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->node.shape.x != x || m_prev->node.shape.y != y);
m_prev->node.shape.left= point;
point->node.shape.right= m_prev;
}
else
m_first= point;
m_prev= point;
m_prev_hook= hook;
return 0;
}
void Gcalc_shape_transporter::int_complete()
{
GCALC_DBUG_ASSERT(m_shape_started == 1 || m_shape_started == 3);
if (!m_first)
return;
/* simple point */
if (m_first == m_prev)
{
m_first->node.shape.right= m_first->node.shape.left= NULL;
return;
}
/* line */
if (m_shape_started == 1)
{
m_first->node.shape.right= NULL;
m_prev->node.shape.left= m_prev->node.shape.right;
m_prev->node.shape.right= NULL;
return;
}
/* polygon */
if (cmp_point_info(m_first, m_prev) == 0)
{
/* Coinciding points, remove the last one from the list */
m_prev->node.shape.right->node.shape.left= m_first;
m_first->node.shape.right= m_prev->node.shape.right;
m_heap->free_point_info(m_prev, m_prev_hook);
}
else
{
GCALC_DBUG_ASSERT(m_prev->node.shape.x != m_first->node.shape.x || m_prev->node.shape.y != m_first->node.shape.y);
m_first->node.shape.right= m_prev;
m_prev->node.shape.left= m_first;
}
}
inline void calc_dx_dy(Gcalc_scan_iterator::point *p)
{
gcalc_sub_coord1(p->dx, p->next_pi->node.shape.ix, p->pi->node.shape.ix);
gcalc_sub_coord1(p->dy, p->next_pi->node.shape.iy, p->pi->node.shape.iy);
if (GCALC_SIGN(p->dx[0]))
{
p->l_border= &p->next_pi->node.shape.ix;
p->r_border= &p->pi->node.shape.ix;
}
else
{
p->r_border= &p->next_pi->node.shape.ix;
p->l_border= &p->pi->node.shape.ix;
}
}
Gcalc_scan_iterator::Gcalc_scan_iterator(size_t blk_size) :
Gcalc_dyn_list(blk_size, sizeof(point) > sizeof(intersection_info) ?
sizeof(point) :
sizeof(intersection_info))
{
state.slice= NULL;
m_bottom_points= NULL;
m_bottom_hook= &m_bottom_points;
}
void Gcalc_scan_iterator::init(Gcalc_heap *points)
{
GCALC_DBUG_ASSERT(points->ready());
GCALC_DBUG_ASSERT(!state.slice);
if (!(m_cur_pi= points->get_first()))
return;
m_heap= points;
state.event_position_hook= &state.slice;
state.event_end= NULL;
#ifndef GCALC_DBUG_OFF
m_cur_thread= 0;
#endif /*GCALC_DBUG_OFF*/
GCALC_SET_TERMINATED(killed, 0);
}
void Gcalc_scan_iterator::reset()
{
state.slice= NULL;
m_bottom_points= NULL;
m_bottom_hook= &m_bottom_points;
Gcalc_dyn_list::reset();
}
int Gcalc_scan_iterator::point::cmp_dx_dy(const Gcalc_coord1 dx_a,
const Gcalc_coord1 dy_a,
const Gcalc_coord1 dx_b,
const Gcalc_coord1 dy_b)
{
Gcalc_coord2 dx_a_dy_b;
Gcalc_coord2 dy_a_dx_b;
gcalc_mul_coord1(dx_a_dy_b, dx_a, dy_b);
gcalc_mul_coord1(dy_a_dx_b, dy_a, dx_b);
return gcalc_cmp_coord(dx_a_dy_b, dy_a_dx_b, GCALC_COORD_BASE2);
}
int Gcalc_scan_iterator::point::cmp_dx_dy(const Gcalc_heap::Info *p1,
const Gcalc_heap::Info *p2,
const Gcalc_heap::Info *p3,
const Gcalc_heap::Info *p4)
{
Gcalc_coord1 dx_a, dy_a, dx_b, dy_b;
gcalc_sub_coord1(dx_a, p2->node.shape.ix, p1->node.shape.ix);
gcalc_sub_coord1(dy_a, p2->node.shape.iy, p1->node.shape.iy);
gcalc_sub_coord1(dx_b, p4->node.shape.ix, p3->node.shape.ix);
gcalc_sub_coord1(dy_b, p4->node.shape.iy, p3->node.shape.iy);
return cmp_dx_dy(dx_a, dy_a, dx_b, dy_b);
}
int Gcalc_scan_iterator::point::cmp_dx_dy(const point *p) const
{
GCALC_DBUG_ASSERT(!is_bottom());
return cmp_dx_dy(dx, dy, p->dx, p->dy);
}
#ifdef GCALC_CHECK_WITH_FLOAT
void Gcalc_scan_iterator::point::calc_x(long double *x, long double y,
long double ix) const
{
long double ddy= gcalc_get_double(dy, GCALC_COORD_BASE);
if (fabsl(ddy) < (long double) 1e-20)
{
*x= ix;
}
else
*x= (ddy * (long double) pi->node.shape.x + gcalc_get_double(dx, GCALC_COORD_BASE) *
(y - pi->node.shape.y)) / ddy;
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
static int compare_events(const void *e0, const void *e1)
{
const Gcalc_scan_iterator::point *p0= (const Gcalc_scan_iterator::point *)e0;
const Gcalc_scan_iterator::point *p1= (const Gcalc_scan_iterator::point *)e1;
return p0->cmp_dx_dy(p1) > 0;
}
int Gcalc_scan_iterator::arrange_event(int do_sorting, int n_intersections)
{
int ev_counter;
point *sp;
point **sp_hook;
ev_counter= 0;
*m_bottom_hook= NULL;
for (sp= m_bottom_points; sp; sp= sp->get_next())
sp->ev_next= sp->get_next();
for (sp= state.slice, sp_hook= &state.slice;
sp; sp_hook= sp->next_ptr(), sp= sp->get_next())
{
if (sp->event)
{
state.event_position_hook= sp_hook;
break;
}
}
for (sp= *(sp_hook= state.event_position_hook);
sp && sp->event; sp_hook= sp->next_ptr(), sp= sp->get_next())
{
ev_counter++;
if (sp->get_next() && sp->get_next()->event)
sp->ev_next= sp->get_next();
else
sp->ev_next= m_bottom_points;
}
#ifndef GCALC_DBUG_OFF
{
point *cur_p= sp;
for (; cur_p; cur_p= cur_p->get_next())
GCALC_DBUG_ASSERT(!cur_p->event);
}
#endif /*GCALC_DBUG_OFF*/
state.event_end= sp;
if (ev_counter == 2 && n_intersections == 1)
{
/* If we had only intersection, just swap the two points. */
sp= *state.event_position_hook;
*state.event_position_hook= sp->get_next();
sp->next= (*state.event_position_hook)->next;
(*state.event_position_hook)->next= sp;
/* The list of the events should be restored. */
(*state.event_position_hook)->ev_next= sp;
sp->ev_next= m_bottom_points;
}
else if (ev_counter == 2 && get_events()->event == scev_two_threads)
{
/* Do nothing. */
}
else if (ev_counter > 1 && do_sorting)
{
point *cur_p;
*sp_hook= NULL;
sp= (point *) sort_list(compare_events, *state.event_position_hook,
ev_counter);
/* Find last item in the list, it's changed after the sorting. */
for (cur_p= sp->get_next(); cur_p->get_next();
cur_p= cur_p->get_next())
{}
cur_p->next= state.event_end;
*state.event_position_hook= sp;
/* The list of the events should be restored. */
for (; sp && sp->event; sp= sp->get_next())
{
if (sp->get_next() && sp->get_next()->event)
sp->ev_next= sp->get_next();
else
sp->ev_next= m_bottom_points;
}
}
#ifndef GCALC_DBUG_OFF
{
const event_point *ev= get_events();
for (; ev && ev->get_next(); ev= ev->get_next())
{
if (ev->is_bottom() || ev->get_next()->is_bottom())
break;
GCALC_DBUG_ASSERT(ev->cmp_dx_dy(ev->get_next()) <= 0);
}
}
#endif /*GCALC_DBUG_OFF*/
return 0;
}
int Gcalc_heap::Info::equal_pi(const Info *pi) const
{
if (type == nt_intersection)
return node.intersection.equal;
if (pi->type == nt_eq_node)
return 1;
if (type == nt_eq_node || pi->type == nt_intersection)
return 0;
return cmp_point_info(this, pi) == 0;
}
int Gcalc_scan_iterator::step()
{
int result= 0;
int do_sorting= 0;
int n_intersections= 0;
point *sp;
GCALC_DBUG_ENTER("Gcalc_scan_iterator::step");
GCALC_DBUG_ASSERT(more_points());
if (GCALC_TERMINATED(killed))
GCALC_DBUG_RETURN(0xFFFF);
/* Clear the old event marks. */
if (m_bottom_points)
{
free_list((Gcalc_dyn_list::Item **) &m_bottom_points,
(Gcalc_dyn_list::Item **) m_bottom_hook);
m_bottom_points= NULL;
m_bottom_hook= &m_bottom_points;
}
for (sp= *state.event_position_hook;
sp != state.event_end; sp= sp->get_next())
sp->event= scev_none;
//#ifndef GCALC_DBUG_OFF
state.event_position_hook= NULL;
state.pi= NULL;
//#endif /*GCALC_DBUG_OFF*/
do
{
#ifndef GCALC_DBUG_OFF
if (m_cur_pi->type == Gcalc_heap::nt_intersection &&
m_cur_pi->get_next()->type == Gcalc_heap::nt_intersection &&
m_cur_pi->node.intersection.equal)
GCALC_DBUG_ASSERT(cmp_intersections(m_cur_pi, m_cur_pi->get_next()) == 0);
#endif /*GCALC_DBUG_OFF*/
GCALC_DBUG_CHECK_COUNTER();
GCALC_DBUG_PRINT_SLICE("step:", state.slice);
GCALC_DBUG_PRINT_PI(m_cur_pi);
if (m_cur_pi->type == Gcalc_heap::nt_shape_node)
{
if (m_cur_pi->is_top())
{
result= insert_top_node();
if (!m_cur_pi->is_bottom())
do_sorting++;
}
else if (m_cur_pi->is_bottom())
remove_bottom_node();
else
{
do_sorting++;
result= node_scan();
}
if (result)
GCALC_DBUG_RETURN(result);
state.pi= m_cur_pi;
}
else if (m_cur_pi->type == Gcalc_heap::nt_eq_node)
{
do_sorting++;
eq_scan();
}
else
{
/* nt_intersection */
do_sorting++;
n_intersections++;
intersection_scan();
if (!state.pi || state.pi->type == Gcalc_heap::nt_intersection)
state.pi= m_cur_pi;
}
m_cur_pi= m_cur_pi->get_next();
} while (m_cur_pi && state.pi->equal_pi(m_cur_pi));
GCALC_DBUG_RETURN(arrange_event(do_sorting, n_intersections));
}
static int node_on_right(const Gcalc_heap::Info *node,
const Gcalc_heap::Info *edge_a, const Gcalc_heap::Info *edge_b)
{
Gcalc_coord1 a_x, a_y;
Gcalc_coord1 b_x, b_y;
Gcalc_coord2 ax_by, ay_bx;
int result;
gcalc_sub_coord1(a_x, node->node.shape.ix, edge_a->node.shape.ix);
gcalc_sub_coord1(a_y, node->node.shape.iy, edge_a->node.shape.iy);
gcalc_sub_coord1(b_x, edge_b->node.shape.ix, edge_a->node.shape.ix);
gcalc_sub_coord1(b_y, edge_b->node.shape.iy, edge_a->node.shape.iy);
gcalc_mul_coord1(ax_by, a_x, b_y);
gcalc_mul_coord1(ay_bx, a_y, b_x);
result= gcalc_cmp_coord(ax_by, ay_bx, GCALC_COORD_BASE2);
#ifdef GCALC_CHECK_WITH_FLOAT
{
long double dx= gcalc_get_double(edge_b->node.shape.ix, GCALC_COORD_BASE) -
gcalc_get_double(edge_a->node.shape.ix, GCALC_COORD_BASE);
long double dy= gcalc_get_double(edge_b->node.shape.iy, GCALC_COORD_BASE) -
gcalc_get_double(edge_a->node.shape.iy, GCALC_COORD_BASE);
long double ax= gcalc_get_double(node->node.shape.ix, GCALC_COORD_BASE) -
gcalc_get_double(edge_a->node.shape.ix, GCALC_COORD_BASE);
long double ay= gcalc_get_double(node->node.shape.iy, GCALC_COORD_BASE) -
gcalc_get_double(edge_a->node.shape.iy, GCALC_COORD_BASE);
long double d= ax * dy - ay * dx;
if (result == 0)
GCALC_DBUG_ASSERT(de_check(d, 0.0));
else if (result < 0)
GCALC_DBUG_ASSERT(de_check(d, 0.0) || d < 0);
else
GCALC_DBUG_ASSERT(de_check(d, 0.0) || d > 0);
}
#endif /*GCALC_CHECK_WITH_FLOAT*/
return result;
}
static int cmp_tops(const Gcalc_heap::Info *top_node,
const Gcalc_heap::Info *edge_a, const Gcalc_heap::Info *edge_b)
{
int cmp_res_a, cmp_res_b;
cmp_res_a= gcalc_cmp_coord1(edge_a->node.shape.ix, top_node->node.shape.ix);
cmp_res_b= gcalc_cmp_coord1(edge_b->node.shape.ix, top_node->node.shape.ix);
if (cmp_res_a <= 0 && cmp_res_b > 0)
return -1;
if (cmp_res_b <= 0 && cmp_res_a > 0)
return 1;
if (cmp_res_a == 0 && cmp_res_b == 0)
return 0;
return node_on_right(edge_a, top_node, edge_b);
}
int Gcalc_scan_iterator::insert_top_node()
{
point *sp= state.slice;
point **prev_hook= &state.slice;
point *sp1= NULL;
point *sp0= new_slice_point();
int cmp_res;
GCALC_DBUG_ENTER("Gcalc_scan_iterator::insert_top_node");
if (!sp0)
GCALC_DBUG_RETURN(1);
sp0->pi= m_cur_pi;
sp0->next_pi= m_cur_pi->node.shape.left;
#ifndef GCALC_DBUG_OFF
sp0->thread= m_cur_thread++;
#endif /*GCALC_DBUG_OFF*/
if (m_cur_pi->node.shape.left)
{
calc_dx_dy(sp0);
if (m_cur_pi->node.shape.right)
{
if (!(sp1= new_slice_point()))
GCALC_DBUG_RETURN(1);
sp1->event= sp0->event= scev_two_threads;
sp1->pi= m_cur_pi;
sp1->next_pi= m_cur_pi->node.shape.right;
#ifndef GCALC_DBUG_OFF
sp1->thread= m_cur_thread++;
#endif /*GCALC_DBUG_OFF*/
calc_dx_dy(sp1);
/* We have two threads so should decide which one will be first */
cmp_res= cmp_tops(m_cur_pi, m_cur_pi->node.shape.left, m_cur_pi->node.shape.right);
if (cmp_res > 0)
{
point *tmp= sp0;
sp0= sp1;
sp1= tmp;
}
else if (cmp_res == 0)
{
/* Exactly same direction of the edges. */
cmp_res= gcalc_cmp_coord1(m_cur_pi->node.shape.left->node.shape.iy, m_cur_pi->node.shape.right->node.shape.iy);
if (cmp_res != 0)
{
if (cmp_res < 0)
{
if (add_eq_node(sp0->next_pi, sp1))
GCALC_DBUG_RETURN(1);
}
else
{
if (add_eq_node(sp1->next_pi, sp0))
GCALC_DBUG_RETURN(1);
}
}
else
{
cmp_res= gcalc_cmp_coord1(m_cur_pi->node.shape.left->node.shape.ix, m_cur_pi->node.shape.right->node.shape.ix);
if (cmp_res != 0)
{
if (cmp_res < 0)
{
if (add_eq_node(sp0->next_pi, sp1))
GCALC_DBUG_RETURN(1);
}
else
{
if (add_eq_node(sp1->next_pi, sp0))
GCALC_DBUG_RETURN(1);
}
}
}
}
}
else
sp0->event= scev_thread;
}
else
sp0->event= scev_single_point;
/* Check if we already have an event - then we'll place the node there */
for (; sp && !sp->event; prev_hook= sp->next_ptr(), sp=sp->get_next())
{}
if (!sp)
{
sp= state.slice;
prev_hook= &state.slice;
/* We need to find the place to insert. */
for (; sp; prev_hook= sp->next_ptr(), sp=sp->get_next())
{
if (sp->event || gcalc_cmp_coord1(*sp->r_border, m_cur_pi->node.shape.ix) < 0)
continue;
cmp_res= node_on_right(m_cur_pi, sp->pi, sp->next_pi);
if (cmp_res == 0)
{
/* The top node lies on the edge. */
/* Nodes of that edge will be handled in other places. */
sp->event= scev_intersection;
}
else if (cmp_res < 0)
break;
}
}
if (sp0->event == scev_single_point)
{
/* Add single point to the bottom list. */
*m_bottom_hook= sp0;
m_bottom_hook= sp0->next_ptr();
state.event_position_hook= prev_hook;
}
else
{
*prev_hook= sp0;
sp0->next= sp;
if (add_events_for_node(sp0))
GCALC_DBUG_RETURN(1);
if (sp0->event == scev_two_threads)
{
*prev_hook= sp1;
sp1->next= sp;
if (add_events_for_node(sp1))
GCALC_DBUG_RETURN(1);
sp0->next= sp1;
*prev_hook= sp0;
}
}
GCALC_DBUG_RETURN(0);
}
void Gcalc_scan_iterator::remove_bottom_node()
{
point *sp= state.slice;
point **sp_hook= &state.slice;
point *first_bottom_point= NULL;
GCALC_DBUG_ENTER("Gcalc_scan_iterator::remove_bottom_node");
for (; sp; sp= sp->get_next())
{
if (sp->next_pi == m_cur_pi)
{
*sp_hook= sp->get_next();
sp->pi= m_cur_pi;
sp->next_pi= NULL;
if (first_bottom_point)
{
first_bottom_point->event= sp->event= scev_two_ends;
break;
}
first_bottom_point= sp;
sp->event= scev_end;
state.event_position_hook= sp_hook;
}
else
sp_hook= sp->next_ptr();
}
GCALC_DBUG_ASSERT(first_bottom_point);
*m_bottom_hook= first_bottom_point;
m_bottom_hook= first_bottom_point->next_ptr();
if (sp)
{
*m_bottom_hook= sp;
m_bottom_hook= sp->next_ptr();
}
GCALC_DBUG_VOID_RETURN;
}
int Gcalc_scan_iterator::add_events_for_node(point *sp_node)
{
point *sp= state.slice;
int cur_pi_r, sp_pi_r;
GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_events_for_node");
/* Scan to the event point. */
for (; sp != sp_node; sp= sp->get_next())
{
GCALC_DBUG_ASSERT(!sp->is_bottom());
GCALC_DBUG_PRINT(("left cut_edge %d", sp->thread));
if (sp->next_pi == sp_node->next_pi ||
gcalc_cmp_coord1(*sp->r_border, *sp_node->l_border) < 0)
continue;
sp_pi_r= node_on_right(sp->next_pi, sp_node->pi, sp_node->next_pi);
if (sp_pi_r < 0)
continue;
cur_pi_r= node_on_right(sp_node->next_pi, sp->pi, sp->next_pi);
if (cur_pi_r > 0)
continue;
if (cur_pi_r == 0 && sp_pi_r == 0)
{
int cmp_res= cmp_point_info(sp->next_pi, sp_node->next_pi);
if (cmp_res > 0)
{
if (add_eq_node(sp_node->next_pi, sp))
GCALC_DBUG_RETURN(1);
}
else if (cmp_res < 0)
{
if (add_eq_node(sp->next_pi, sp_node))
GCALC_DBUG_RETURN(1);
}
continue;
}
if (cur_pi_r == 0)
{
if (add_eq_node(sp_node->next_pi, sp))
GCALC_DBUG_RETURN(1);
continue;
}
else if (sp_pi_r == 0)
{
if (add_eq_node(sp->next_pi, sp_node))
GCALC_DBUG_RETURN(1);
continue;
}
if (sp->event)
{
#ifndef GCALC_DBUG_OFF
cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi);
GCALC_DBUG_ASSERT(cur_pi_r == 0);
#endif /*GCALC_DBUG_OFF*/
continue;
}
cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi);
GCALC_DBUG_ASSERT(cur_pi_r >= 0);
//GCALC_DBUG_ASSERT(cur_pi_r > 0); /* Is it ever violated? */
if (cur_pi_r > 0 && add_intersection(sp, sp_node, m_cur_pi))
GCALC_DBUG_RETURN(1);
}
/* Scan to the end of the slice */
sp= sp->get_next();
for (; sp; sp= sp->get_next())
{
GCALC_DBUG_ASSERT(!sp->is_bottom());
GCALC_DBUG_PRINT(("right cut_edge %d", sp->thread));
if (sp->next_pi == sp_node->next_pi ||
gcalc_cmp_coord1(*sp_node->r_border, *sp->l_border) < 0)
continue;
sp_pi_r= node_on_right(sp->next_pi, sp_node->pi, sp_node->next_pi);
if (sp_pi_r > 0)
continue;
cur_pi_r= node_on_right(sp_node->next_pi, sp->pi, sp->next_pi);
if (cur_pi_r < 0)
continue;
if (cur_pi_r == 0 && sp_pi_r == 0)
{
int cmp_res= cmp_point_info(sp->next_pi, sp_node->next_pi);
if (cmp_res > 0)
{
if (add_eq_node(sp_node->next_pi, sp))
GCALC_DBUG_RETURN(1);
}
else if (cmp_res < 0)
{
if (add_eq_node(sp->next_pi, sp_node))
GCALC_DBUG_RETURN(1);
}
continue;
}
if (cur_pi_r == 0)
{
if (add_eq_node(sp_node->next_pi, sp))
GCALC_DBUG_RETURN(1);
continue;
}
else if (sp_pi_r == 0)
{
if (add_eq_node(sp->next_pi, sp_node))
GCALC_DBUG_RETURN(1);
continue;
}
if (sp->event)
{
#ifndef GCALC_DBUG_OFF
cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi);
GCALC_DBUG_ASSERT(cur_pi_r == 0);
#endif /*GCALC_DBUG_OFF*/
continue;
}
cur_pi_r= node_on_right(sp_node->pi, sp->pi, sp->next_pi);
GCALC_DBUG_ASSERT(cur_pi_r <= 0);
//GCALC_DBUG_ASSERT(cur_pi_r < 0); /* Is it ever violated? */
if (cur_pi_r < 0 && add_intersection(sp_node, sp, m_cur_pi))
GCALC_DBUG_RETURN(1);
}
GCALC_DBUG_RETURN(0);
}
int Gcalc_scan_iterator::node_scan()
{
point *sp= state.slice;
Gcalc_heap::Info *cur_pi= m_cur_pi;
GCALC_DBUG_ENTER("Gcalc_scan_iterator::node_scan");
/* Scan to the event point. */
/* Can be avoided if we add link to the sp to the Info. */
for (; sp->next_pi != cur_pi; sp= sp->get_next())
{}
GCALC_DBUG_PRINT(("node for %d", sp->thread));
/* Handle the point itself. */
sp->pi= cur_pi;
sp->next_pi= cur_pi->node.shape.left;
sp->event= scev_point;
calc_dx_dy(sp);
GCALC_DBUG_RETURN(add_events_for_node(sp));
}
void Gcalc_scan_iterator::eq_scan()
{
point *sp= eq_sp(m_cur_pi);
GCALC_DBUG_ENTER("Gcalc_scan_iterator::eq_scan");
#ifndef GCALC_DBUG_OFF
{
point *cur_p= state.slice;
for (; cur_p && cur_p != sp; cur_p= cur_p->get_next())
{}
GCALC_DBUG_ASSERT(cur_p);
}
#endif /*GCALC_DBUG_OFF*/
if (!sp->event)
{
sp->event= scev_intersection;
sp->ev_pi= m_cur_pi;
}
GCALC_DBUG_VOID_RETURN;
}
void Gcalc_scan_iterator::intersection_scan()
{
intersection_info *ii= i_data(m_cur_pi);
GCALC_DBUG_ENTER("Gcalc_scan_iterator::intersection_scan");
#ifndef GCALC_DBUG_OFF
{
point *sp= state.slice;
for (; sp && sp != ii->edge_a; sp= sp->get_next())
{}
GCALC_DBUG_ASSERT(sp);
for (; sp && sp != ii->edge_b; sp= sp->get_next())
{}
GCALC_DBUG_ASSERT(sp);
}
#endif /*GCALC_DBUG_OFF*/
ii->edge_a->event= ii->edge_b->event= scev_intersection;
ii->edge_a->ev_pi= ii->edge_b->ev_pi= m_cur_pi;
free_item(ii);
m_cur_pi->node.intersection.data= NULL;
GCALC_DBUG_VOID_RETURN;
}
int Gcalc_scan_iterator::add_intersection(point *sp_a, point *sp_b,
Gcalc_heap::Info *pi_from)
{
Gcalc_heap::Info *ii;
intersection_info *i_calc;
int cmp_res;
int skip_next= 0;
GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_intersection");
if (!(i_calc= new_intersection_info(sp_a, sp_b)) ||
!(ii= new_intersection(m_heap, i_calc)))
GCALC_DBUG_RETURN(1);
ii->node.intersection.equal= 0;
for (;
pi_from->get_next() != sp_a->next_pi &&
pi_from->get_next() != sp_b->next_pi;
pi_from= pi_from->get_next())
{
Gcalc_heap::Info *cur= pi_from->get_next();
if (skip_next)
{
if (cur->type == Gcalc_heap::nt_intersection)
skip_next= cur->node.intersection.equal;
else
skip_next= 0;
continue;
}
if (cur->type == Gcalc_heap::nt_intersection)
{
cmp_res= cmp_intersections(cur, ii);
skip_next= cur->node.intersection.equal;
}
else if (cur->type == Gcalc_heap::nt_eq_node)
continue;
else
cmp_res= cmp_node_isc(cur, ii);
if (cmp_res == 0)
{
ii->node.intersection.equal= 1;
break;
}
else if (cmp_res > 0)
break;
}
/* Intersection inserted before the equal point. */
ii->next= pi_from->get_next();
pi_from->next= ii;
GCALC_DBUG_RETURN(0);
}
int Gcalc_scan_iterator::add_eq_node(Gcalc_heap::Info *node, point *sp)
{
Gcalc_heap::Info *en;
GCALC_DBUG_ENTER("Gcalc_scan_iterator::add_intersection");
en= new_eq_point(m_heap, node, sp);
if (!en)
GCALC_DBUG_RETURN(1);
/* eq_node iserted after teh equal point. */
en->next= node->get_next();
node->next= en;
GCALC_DBUG_RETURN(0);
}
void calc_t(Gcalc_coord2 t_a, Gcalc_coord2 t_b,
Gcalc_coord1 dxa, Gcalc_coord1 dxb,
const Gcalc_heap::Info *p1, const Gcalc_heap::Info *p2,
const Gcalc_heap::Info *p3, const Gcalc_heap::Info *p4)
{
Gcalc_coord1 a2_a1x, a2_a1y;
Gcalc_coord2 x1y2, x2y1;
Gcalc_coord1 dya, dyb;
gcalc_sub_coord1(a2_a1x, p3->node.shape.ix, p1->node.shape.ix);
gcalc_sub_coord1(a2_a1y, p3->node.shape.iy, p1->node.shape.iy);
gcalc_sub_coord1(dxa, p2->node.shape.ix, p1->node.shape.ix);
gcalc_sub_coord1(dya, p2->node.shape.iy, p1->node.shape.iy);
gcalc_sub_coord1(dxb, p4->node.shape.ix, p3->node.shape.ix);
gcalc_sub_coord1(dyb, p4->node.shape.iy, p3->node.shape.iy);
gcalc_mul_coord1(x1y2, dxa, dyb);
gcalc_mul_coord1(x2y1, dya, dxb);
gcalc_sub_coord(t_b, GCALC_COORD_BASE2, x1y2, x2y1);
gcalc_mul_coord1(x1y2, a2_a1x, dyb);
gcalc_mul_coord1(x2y1, a2_a1y, dxb);
gcalc_sub_coord(t_a, GCALC_COORD_BASE2, x1y2, x2y1);
}
double Gcalc_scan_iterator::get_y() const
{
if (state.pi->type == Gcalc_heap::nt_intersection)
{
Gcalc_coord1 dxa, dya;
Gcalc_coord2 t_a, t_b;
Gcalc_coord3 a_tb, b_ta, y_exp;
calc_t(t_a, t_b, dxa, dya,
state.pi->node.intersection.p1, state.pi->node.intersection.p2, state.pi->node.intersection.p3, state.pi->node.intersection.p4);
gcalc_mul_coord(a_tb, GCALC_COORD_BASE3,
t_b, GCALC_COORD_BASE2, state.pi->node.intersection.p1->node.shape.iy, GCALC_COORD_BASE);
gcalc_mul_coord(b_ta, GCALC_COORD_BASE3,
t_a, GCALC_COORD_BASE2, dya, GCALC_COORD_BASE);
gcalc_add_coord(y_exp, GCALC_COORD_BASE3, a_tb, b_ta);
return (get_pure_double(y_exp, GCALC_COORD_BASE3) /
get_pure_double(t_b, GCALC_COORD_BASE2)) / m_heap->coord_extent;
}
else
return state.pi->node.shape.y;
}
double Gcalc_scan_iterator::get_event_x() const
{
if (state.pi->type == Gcalc_heap::nt_intersection)
{
Gcalc_coord1 dxa, dya;
Gcalc_coord2 t_a, t_b;
Gcalc_coord3 a_tb, b_ta, x_exp;
calc_t(t_a, t_b, dxa, dya,
state.pi->node.intersection.p1, state.pi->node.intersection.p2, state.pi->node.intersection.p3, state.pi->node.intersection.p4);
gcalc_mul_coord(a_tb, GCALC_COORD_BASE3,
t_b, GCALC_COORD_BASE2, state.pi->node.intersection.p1->node.shape.ix, GCALC_COORD_BASE);
gcalc_mul_coord(b_ta, GCALC_COORD_BASE3,
t_a, GCALC_COORD_BASE2, dxa, GCALC_COORD_BASE);
gcalc_add_coord(x_exp, GCALC_COORD_BASE3, a_tb, b_ta);
return (get_pure_double(x_exp, GCALC_COORD_BASE3) /
get_pure_double(t_b, GCALC_COORD_BASE2)) / m_heap->coord_extent;
}
else
return state.pi->node.shape.x;
}
double Gcalc_scan_iterator::get_h() const
{
double cur_y= get_y();
double next_y;
if (state.pi->type == Gcalc_heap::nt_intersection)
{
double x;
state.pi->calc_xy(&x, &next_y);
}
else
next_y= state.pi->next ? state.pi->get_next()->node.shape.y : 0.0;
return next_y - cur_y;
}
double Gcalc_scan_iterator::get_sp_x(const point *sp) const
{
double dy;
if (sp->event & (scev_end | scev_two_ends | scev_point))
return sp->pi->node.shape.x;
dy= sp->next_pi->node.shape.y - sp->pi->node.shape.y;
if (fabs(dy) < 1e-12)
return sp->pi->node.shape.x;
return sp->pi->node.shape.x + (sp->next_pi->node.shape.x - sp->pi->node.shape.x) * dy;
}
double Gcalc_scan_iterator::get_pure_double(const Gcalc_internal_coord *d,
int d_len)
{
int n= 1;
long double res= (long double) FIRST_DIGIT(d[0]);
do
{
res*= (long double) GCALC_DIG_BASE;
res+= (long double) d[n];
} while(++n < d_len);
if (GCALC_SIGN(d[0]))
res*= -1.0;
return res;
}
#endif /* HAVE_SPATIAL */
Reference in a new issue View git blame Copy permalink