44 void heapsort(
int *sid,
int *
id,
const double*
const x,
int N )
46 unsigned int n = N, i = n / 2, parent, child;
66 if ( child + 1 < n && x[
id[sid[child + 1]]] > x[
id[sid[child]]] )
70 if ( x[
id[sid[child]]] > x[
id[tx]] )
72 sid[parent] = sid[child];
74 child = parent * 2 + 1;
88 unsigned int n = N, i = n / 2, parent, child;
102 if ( n == 0 )
return;
112 if ( child + 1 < n && heap[child + 1] > heap[child] )
116 if ( heap[child] > t )
118 heap[parent] = heap[child];
119 x[parent] = x[child];
121 child = parent * 2 + 1;
141 double x3,
double y3,
double x4,
double y4 )
150 return (
cross_product( x1, y1, x2, y2, x3, y3 ) *
cross_product( x1, y1, x2, y2, x4, y4 ) < 0
151 &&
cross_product( x3, y3, x4, y4, x1, y1 ) *
cross_product( x3, y3, x4, y4, x2, y2 ) < 0 );
161 double x3,
double y3,
double x4,
double y4,
double xs2,
double ys2,
162 double *x,
double *y )
164 double cp1, cp2, cp3, cp4;
171 if ( cp1 == 0 && cp2 == 0 && cp3 == 0 && cp4 == 0 )
174 std::cout <<
"coolineaire..." << std::endl;
180 if ( cp1 == 0 && cp3 == 0 )
183 std::cout <<
"cp1 = cp3 = 0 => ignoring..." << std::endl;
189 if ( cp1 == 0 && cp4 == 0 )
192 std::cout <<
"cp1 = cp4 = 0 => ignoring..." << std::endl;
198 if ( cp2 == 0 && cp3 == 0 )
201 std::cout <<
"cp2 = cp3 = 0 => ignoring..." << std::endl;
207 if ( cp1 == 0 || cp3 == 0 )
210 std::cout <<
"skip..." << std::endl;
216 if ( cp4 == 0 && cp1 * cp1 < 0 )
229 if ( cp2 == 0 && cp3 * cp4 < 0 )
242 if ( cp2 == 0 && cp4 == 0 )
268 toDist = distance[0];
271 toDist =
max( toDist, distance[1] );
274 toDist =
max( toDist, distance[2] );
277 toDist =
max( toDist, distance[3] );
279 for ( i = 0; i < 4; i++ )
284 ratio = toDist / distance[i];
286 nx[i] = cx + dx * ratio;
287 ny[i] = cy + dy * ratio;
312 double x3,
double y3,
double x4,
double y4,
313 double *x,
double *y )
319 if ( cp1 * cp2 <= 0 )
333 double x3,
double y3,
double x4,
double y4,
334 double *x,
double *y )
336 double cp1, cp2, cp3, cp4;
354 double x3,
double y3,
double x4,
double y4,
355 double *x,
double *y )
358 double a1, a2, b1, b2, c1, c2;
363 c1 = x2 * y1 - x1 * y2;
367 c2 = x4 * y3 - x3 * y4;
369 if (( denom = a1 * b2 - a2 * b1 ) == 0 )
375 *x = ( b1 * c2 - b2 * c1 ) / denom;
376 *y = ( a2 * c1 - a1 * c2 ) / denom;
394 int convexHullId(
int *
id,
const double*
const x,
const double*
const y,
int n,
int *&cHull )
399 for ( i = 0; i < n; i++ )
405 if ( n <= 3 )
return n;
407 int* stack =
new int[n];
408 double* tan =
new double [n];
419 while ( ref < n &&
vabs( y[
id[cHull[ref]]] - y[
id[cHull[0]]] ) <
EPSILON ) ref++;
424 for ( i = ref; i < n; i++ )
426 if (
vabs( y[
id[cHull[i]]] - y[
id[cHull[0]]] ) <
EPSILON )
429 tan[i] = ( x[
id[cHull[0]]] - x[
id[cHull[i]]] ) / ( y[
id[cHull[i]]] - y[
id[cHull[0]]] );
433 heapsort2( cHull + ref, tan + ref, n - ref );
443 stack[1] = cHull[ref-1];
449 for ( i = ref; i < n; i++ )
451 result =
cross_product( x[
id[stack[second]]], y[
id[stack[second]]],
452 x[
id[stack[top]]], y[
id[stack[top]]], x[
id[cHull[i]]], y[
id[cHull[i]]] );
456 if (
dist_euc2d_sq( x[
id[stack[second]]], y[
id[stack[second]]], x[
id[cHull[i]]], y[
id[cHull[i]]] )
457 >
dist_euc2d_sq( x[
id[stack[second]]], y[
id[stack[second]]], x[
id[stack[top]]], y[
id[stack[top]]] ) )
459 stack[top] = cHull[i];
462 else if ( result > 0 )
466 stack[top] = cHull[i];
470 while ( result < 0 && top > 1 )
475 y[
id[stack[second]]], x[
id[stack[top]]],
476 y[
id[stack[top]]], x[
id[cHull[i]]], y[
id[cHull[i]]] );
480 stack[top] = cHull[i];
484 for ( i = 0; i <= top; i++ )
503 int *pts =
new int[nbPoints];
504 for ( i = 0; i < nbPoints; i++ )
510 if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] )
512 else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] )
514 else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
516 else if ( pts[cHull[0]] > pts[cHull[1]] && pts[cHull[1]] < pts[cHull[2]] && pts[cHull[2]] > pts[cHull[0]] )
518 else if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] > pts[cHull[0]] )
520 else if ( pts[cHull[0]] < pts[cHull[1]] && pts[cHull[1]] > pts[cHull[2]] && pts[cHull[2]] < pts[cHull[0]] )
524 std::cout <<
"Warning wrong cHull -> geometry: " << nbPoints << std::endl;
525 for ( i = 0; i < nbPoints; i++ )
527 std::cout << x[i] <<
";" << y[i] << std::endl;
529 std::cout <<
"hull : " << cHullSize << std::endl;
530 for ( i = 0; i < cHullSize; i++ )
532 std::cout << pts[cHull[i]] <<
" ";
534 std::cout << std::endl;
544 for ( i = 0, j = nbPoints - 1; i <= j; i++, j-- )
571 for ( i = 0, j = npol - 1; i < npol; j = i++ )
573 if (((( yp[i] <= y ) && ( y < yp[j] ) ) ||
574 (( yp[j] <= y ) && ( y < yp[i] ) ) )
575 && ( x < ( xp[j] - xp[i] ) *( y - yp[i] ) / ( yp[j] - yp[i] ) + xp[i] ) )
584 void toSVGPath(
int nbPoints,
int geomType,
double *x,
double *y,
585 int dpi,
double scale,
int xmin,
int ymax,
586 char *layername,
char *objectID,
593 out <<
" <path style=\"fill:none;fill-opacity:1;fill-rule:evenodd;stroke:#000000;stroke-width:1;stroke-linecap:round;stroke-linejoin:round;stroke-opacity:1\" d=\"M " <<
convert2pt( x[0], scale, dpi ) - xmin <<
"," << ymax -
convert2pt( y[0], scale, dpi ) <<
" ";
594 for ( i = 1; i < nbPoints; i++ )
596 out <<
"L " <<
convert2pt( x[i], scale, dpi ) - xmin <<
", " << ymax -
convert2pt( y[i], scale, dpi ) <<
" ";
599 if ( geomType == GEOS_POLYGON )
605 out <<
"id=\"" << layername <<
"-" << objectID <<
"\" ";
606 out <<
"inkscape:label=\"#path-" << layername <<
"-" << objectID <<
"\"/>\n";
610 int cx =
convert2pt( x[0], scale, dpi ) - xmin;
611 int cy = ymax -
convert2pt( y[0], scale, dpi );
613 out <<
" sodipodi:type=\"arc\" ";
614 out <<
" style=\"opacity:1;fill:#bcbcbc;fill-opacity:l;stroke:#000000;stroke-opacity:1;stroke-width:0.5;stroke-linejoin:miter;stroke-dasharray:none;display:inline\"";
615 out <<
" id=\"" << layername <<
"-" << objectID <<
"\" ";
616 out <<
" sodipodi:cx=\"" << cx <<
"\" ";
617 out <<
" sodipodi:cy=\"" << cy <<
"\" ";
618 out <<
" sodipodi:rx=\"1\" ";
619 out <<
" sodipodi:ry=\"1\" ";
620 out <<
" d=\"M " << cx + 1 <<
" " << cy <<
" A 1 1 0 1 1 "
621 << cx - 1 <<
"," << cy <<
" A 1 1 0 1 1 "
622 << cx + 1 <<
" " << cy <<
" z\" ";
623 out <<
" inkscape:label=\"#path-" << layername <<
"-" << objectID <<
"\" />\n";
630 double x1,
double y1,
double x2,
double y2,
631 double& xRes,
double& yRes )
636 double A = dx * dx + dy * dy;
637 double B = 2 * ( dx * ( x1 - cx ) + dy * ( y1 - cy ) );
638 double C = ( x1 - cx ) * ( x1 - cx ) + ( y1 - cy ) * ( y1 - cy ) - radius * radius;
640 double det = B * B - 4 * A * C;
641 if ( A <= 0.0000001 || det < 0 )
648 double t = -B / ( 2 * A );
657 double t = ( -B + sqrt( det ) ) / ( 2 * A );
int reorderPolygon(int nbPoints, double *x, double *y)
void convert2pt(int *x, double scale, int dpi)
bool computeLineIntersection(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y)
void heapsort(int *sid, int *id, const double *const x, int N)
double dist_euc2d_sq(double x1, double y1, double x2, double y2)
bool computeSegIntersection(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y)
bool isPointInPolygon(int npol, double *xp, double *yp, double x, double y)
bool computeSegIntersectionExt(double x1, double y1, double x2, double y2, double xs1, double ys1, double x3, double y3, double x4, double y4, double xs2, double ys2, double *x, double *y)
double dist_euc2d(double x1, double y1, double x2, double y2)
void heapsort2(int *x, double *heap, int N)
double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
bool computeLineSegIntersection(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y)
int convexHullId(int *id, const double *const x, const double *const y, int n, int *&cHull)
bool isSegIntersects(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4)
void findLineCircleIntersection(double cx, double cy, double radius, double x1, double y1, double x2, double y2, double &xRes, double &yRes)