33 #if defined(_VERBOSE_) || (_DEBUG_) || (_DEBUG_FULL_)
129 memcpy(
x, ps.
x,
sizeof(
double )*
nbPoints );
130 memcpy(
y, ps.
y,
sizeof(
double )* nbPoints );
162 bool needClose =
false;
168 GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt,
nbPoints + ( needClose ? 1 : 0 ), 2 );
169 for (
int i = 0; i <
nbPoints; ++i )
171 GEOSCoordSeq_setX_r( geosctxt, coord, i,
x[i] );
172 GEOSCoordSeq_setY_r( geosctxt, coord, i,
y[i] );
178 GEOSCoordSeq_setX_r( geosctxt, coord, nbPoints,
x[0] );
179 GEOSCoordSeq_setY_r( geosctxt, coord, nbPoints,
y[0] );
185 mGeos = GEOSGeom_createPolygon_r( geosctxt, GEOSGeom_createLinearRing_r( geosctxt, coord ), 0, 0 );
188 case GEOS_LINESTRING:
189 mGeos = GEOSGeom_createLineString_r( geosctxt, coord );
193 mGeos = GEOSGeom_createPoint_r( geosctxt, coord );
205 if ( !mPreparedGeom )
209 return mPreparedGeom;
216 GEOSGeom_destroy_r( geosctxt,
mGeos );
217 GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
229 GEOSGeom_destroy_r( geosctxt,
mGeos );
232 GEOSPreparedGeom_destroy_r( geosctxt, mPreparedGeom );
255 newShape->
type = GEOS_POLYGON;
259 newShape->
x =
new double[newShape->
nbPoints];
260 newShape->
y =
new double[newShape->
nbPoints];
262 std::cout <<
"New shape: ";
265 for ( j = 0, i = imin; i != ( imax + 1 ) %
nbPoints; i = ( i + 1 ) %
nbPoints, j++ )
267 newShape->
x[j] =
x[i];
268 newShape->
y[j] =
y[i];
270 std::cout <<
x[i] <<
";" <<
y[i] << std::endl;
277 newShape->
x[j] = fptx;
278 newShape->
y[j] = fpty;
280 std::cout << fptx <<
";" << fpty << std::endl;
285 std::cout <<
"J = " << j <<
"/" << newShape->
nbPoints << std::endl;
286 std::cout << std::endl;
287 std::cout <<
"This: " <<
this << std::endl;
298 GEOSCoordSequence* seq = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
299 GEOSCoordSeq_setX_r( geosctxt, seq, 0, x );
300 GEOSCoordSeq_setY_r( geosctxt, seq, 0, y );
301 GEOSGeometry* point = GEOSGeom_createPoint_r( geosctxt, seq );
302 bool result = ( GEOSPreparedContains_r( geosctxt,
preparedGeom(), point ) == 1 );
303 GEOSGeom_destroy_r( geosctxt, point );
307 catch ( GEOSException &e )
318 GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 5, 2 );
320 GEOSCoordSeq_setX_r( geosctxt, coord, 0, x );
321 GEOSCoordSeq_setY_r( geosctxt, coord, 0, y );
324 double beta = alpha + (
M_PI / 2 );
325 double dx1 = cos( alpha ) * width;
326 double dy1 = sin( alpha ) * width;
327 double dx2 = cos( beta ) * height;
328 double dy2 = sin( beta ) * height;
329 GEOSCoordSeq_setX_r( geosctxt, coord, 1, x + dx1 );
330 GEOSCoordSeq_setY_r( geosctxt, coord, 1, y + dy1 );
331 GEOSCoordSeq_setX_r( geosctxt, coord, 2, x + dx1 + dx2 );
332 GEOSCoordSeq_setY_r( geosctxt, coord, 2, y + dy1 + dy2 );
333 GEOSCoordSeq_setX_r( geosctxt, coord, 3, x + dx2 );
334 GEOSCoordSeq_setY_r( geosctxt, coord, 3, y + dy2 );
338 GEOSCoordSeq_setX_r( geosctxt, coord, 1, x + width );
339 GEOSCoordSeq_setY_r( geosctxt, coord, 1, y );
340 GEOSCoordSeq_setX_r( geosctxt, coord, 2, x + width );
341 GEOSCoordSeq_setY_r( geosctxt, coord, 2, y + height );
342 GEOSCoordSeq_setX_r( geosctxt, coord, 3, x );
343 GEOSCoordSeq_setY_r( geosctxt, coord, 3, y + height );
346 GEOSCoordSeq_setX_r( geosctxt, coord, 4, x );
347 GEOSCoordSeq_setY_r( geosctxt, coord, 4, y );
351 GEOSGeometry* bboxGeos = GEOSGeom_createLinearRing_r( geosctxt, coord );
352 bool result = ( GEOSPreparedContains_r( geosctxt,
preparedGeom(), bboxGeos ) == 1 );
353 GEOSGeom_destroy_r( geosctxt, bboxGeos );
356 catch ( GEOSException &e )
368 std::cout <<
"splitPolygons: " << uid << std::endl;
405 double labelArea = xrm * yrm;
409 while ( shapes_toProcess.
size() > 0 )
412 std::cout <<
"Shape popping()" << std::endl;
422 std::cout <<
"nbp: " << nbp << std::endl;
423 std::cout <<
" PtSet: ";
426 for ( i = 0; i < nbp; i++ )
430 std::cout << x[i] <<
";" << y[i] << std::endl;
435 std::cout << std::endl;
441 cHull = shape->
cHull;
446 std::cout <<
" CHull: ";
449 std::cout << cHull[i] <<
" ";
451 std::cout << std::endl;
461 ihn = ( ihs + 1 ) % cHullSize;
464 ipn = ( ips + 1 ) % nbp;
465 if ( ipn != cHull[ihn] )
470 for ( i = ips; i != cHull[ihn]; i = ( i + 1 ) % nbp )
473 x[cHull[ihn]], y[cHull[ihn]],
483 std::cout <<
"Deeper POint: " << pt <<
" between " << ips <<
" and " << cHull[ihn] << std::endl;
488 base =
dist_euc2d( x[cHull[ihs]], y[cHull[ihs]],
489 x[cHull[ihn]], y[cHull[ihn]] );
497 s = ( base + b + c ) / 2;
498 area = s * ( s - base ) * ( s - b ) * ( s - c );
503 if ( area - bestArea >
EPSILON )
520 bestArea = sqrt( bestArea );
521 double cx, cy, dx, dy, ex, ey, fx, fy, seg_length, ptx = 0, pty = 0, fptx = 0, fpty = 0;
522 int ps = -1, pe = -1, fps = -1, fpe = -1;
523 if ( retainedPt >= 0 && bestArea > labelArea )
526 std::cout <<
"Trou: " << retainedPt << std::endl;
527 std::cout <<
"Hole: " << cHull[holeS] <<
" -> " << cHull[holeE] << std::endl;
528 std::cout <<
"iterate from " << ( cHull[holeE] + 1 ) % nbp <<
" to " << ( cHull[holeS] - 1 + nbp ) % nbp << std::endl;
536 for ( i = ( cHull[holeE] + 1 ) % nbp; i != ( cHull[holeS] - 1 + nbp ) % nbp; i = j )
543 seg_length =
dist_euc2d( x[i], y[i], x[j], y[j] );
544 cx = ( x[i] + x[j] ) / 2.0;
545 cy = ( y[i] + y[j] ) / 2.0;
554 std::cout <<
"D: " << dx <<
" " << dy << std::endl;
555 std::cout <<
"seg_length: " << seg_length << std::endl;
557 if ( seg_length <
EPSILON || qAbs(( b =
cross_product( ex, ey, fx, fy, x[retainedPt], y[retainedPt] ) / ( seg_length ) ) ) > ( seg_length / 2 ) )
559 if (( ex =
dist_euc2d_sq( x[i], y[i], x[retainedPt], y[retainedPt] ) ) < ( ey =
dist_euc2d_sq( x[j], y[j], x[retainedPt], y[retainedPt] ) ) )
574 b =
cross_product( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt] ) / seg_length;
579 if ( !
computeLineIntersection( x[i], y[i], x[j], y[j], x[retainedPt], y[retainedPt], x[retainedPt] - dx, y[retainedPt] + dy, &ptx, &pty ) )
581 std::cout <<
"Oups ... il devrait par tomber la..." << std::endl;
584 std::cout <<
"intersection : " << x[i] <<
" " << y[i] <<
" " << x[j] <<
" " << y[j] <<
" " << x[retainedPt] <<
" " << y[retainedPt] <<
" " << x[retainedPt] - dx <<
" " << y[retainedPt] + dy << std::endl;
585 std::cout <<
" => " << ptx <<
";" << pty << std::endl;
586 std::cout <<
" cxy> " << cx <<
";" << cy << std::endl;
587 std::cout <<
" dxy> " << dx <<
";" << dy << std::endl;
592 double pointX, pointY;
604 for ( k = cHull[holeS]; k != cHull[holeE]; k = ( k + 1 ) % nbp )
608 if (
isSegIntersects( x[retainedPt], y[retainedPt], pointX, pointY, x[k], y[k], x[l], y[l] ) )
617 if ( isValid && b < c )
629 std::cout <<
" cut from " << retainedPt <<
" to ";
631 std::cout <<
"point " << fps << std::endl;
634 std::cout <<
"new point (" << fptx <<
";" << fpty <<
" between " << fps <<
" and " << fpe << std::endl;
639 int imin = retainedPt;
640 int imax = ((( fps < retainedPt && fpe < retainedPt ) || ( fps > retainedPt && fpe > retainedPt ) ) ? qMin( fps, fpe ) : qMax( fps, fpe ) );
642 int nbPtSh1, nbPtSh2;
644 nbPtSh1 = imax - imin + 1 + ( fpe != fps );
646 nbPtSh1 = imax + nbp - imin + 1 + ( fpe != fps );
648 if (( imax == fps ? fpe : fps ) < imin )
649 nbPtSh2 = imin - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
651 nbPtSh2 = imin + nbp - ( imax == fps ? fpe : fps ) + 1 + ( fpe != fps );
655 std::cout <<
"imin: " << imin <<
" imax:" << imax << std::endl;
657 if ( retainedPt == -1 || fps == -1 || fpe == -1 )
660 std::cout << std::endl <<
"Failed to split feature !!! (uid=" << uid <<
")" << std::endl;
666 else if ( imax == imin || nbPtSh1 <= 2 || nbPtSh2 <= 2 || nbPtSh1 == nbp || nbPtSh2 == nbp )
668 shapes_final.
append( shape );
682 std::cout <<
"push back:" << std::endl;
683 for ( i = 0; i < newShape->
nbPoints; i++ )
685 std::cout << newShape->
x[i] <<
";" << newShape->
y[i] << std::endl;
689 shapes_toProcess.
append( newShape );
696 newShape = shape->
extractShape( nbPtSh2, imax, imin, fps, fpe, fptx, fpty );
704 std::cout <<
"push back:" << std::endl;
705 for ( i = 0; i < newShape->
nbPoints; i++ )
707 std::cout << newShape->
x[i] <<
";" << newShape->
y[i] << std::endl;
710 shapes_toProcess.
append( newShape );
719 std::cout <<
"Put shape into shapes_final" << std::endl;
721 shapes_final.
append( shape );
746 double distNearestPoint;
752 double best_area = DBL_MAX;
753 double best_alpha = -1;
755 double best_length = 0;
756 double best_width = 0;
765 std::cout <<
"Compute_chull_bbox" << std::endl;
771 std::cout <<
x[
cHull[i]] <<
";" <<
y[
cHull[i]] << std::endl;
773 if (
x[
cHull[i]] < bbox[0] )
776 if (
x[cHull[i]] > bbox[2] )
777 bbox[2] =
x[cHull[i]];
779 if (
y[cHull[i]] < bbox[1] )
780 bbox[1] =
y[cHull[i]];
782 if (
y[cHull[i]] > bbox[3] )
783 bbox[3] =
y[cHull[i]];
787 dref = bbox[2] - bbox[0];
789 for ( alpha_d = 0; alpha_d < 90; alpha_d++ )
791 alpha = alpha_d *
M_PI / 180.0;
792 d1 = cos( alpha ) * dref;
793 d2 = sin( alpha ) * dref;
814 bb[14] = bb[12] + d2;
815 bb[15] = bb[13] - d1;
818 for ( i = 0; i < 16; i += 4 )
821 alpha_seg = (( i / 4 > 0 ? ( i / 4 ) - 1 : 3 ) ) *
M_PI / 2 + alpha;
833 distNearestPoint = best_cp / dref;
835 d1 = cos( alpha_seg ) * distNearestPoint;
836 d2 = sin( alpha_seg ) * distNearestPoint;
845 width =
cross_product( bb[6], bb[7], bb[4], bb[5], bb[12], bb[13] ) / dref;
846 length =
cross_product( bb[2], bb[3], bb[0], bb[1], bb[8], bb[9] ) / dref;
854 if ( best_area - area >
EPSILON )
860 memcpy( best_bb, bb,
sizeof(
double ) *16 );
868 for ( i = 0; i < 16; i = i + 4 )
871 best_bb[( i+4 ) %16], best_bb[( i+5 ) %16], best_bb[( i+6 ) %16], best_bb[( i+7 ) %16],
872 &finalBb->
x[int ( i/4 )], &finalBb->
y[int ( i/4 )] );
875 finalBb->
alpha = best_alpha;
876 finalBb->
width = best_width;
877 finalBb->
length = best_length;
880 std::cout <<
"FINAL" << std::endl;
881 std::cout <<
"Length : " << best_length << std::endl;
882 std::cout <<
"Width : " << best_width << std::endl;
883 std::cout <<
"Alpha: " << best_alpha <<
" " << best_alpha*180 /
M_PI << std::endl;
884 for ( i = 0; i < 4; i++ )
886 std::cout << finalBb->
x[0] <<
" " << finalBb->
y[0] <<
" ";
888 std::cout << std::endl;
905 GEOSCoordSequence *coord = GEOSCoordSeq_create_r( geosctxt, 1, 2 );
906 GEOSCoordSeq_setX_r( geosctxt, coord, 0, px );
907 GEOSCoordSeq_setY_r( geosctxt, coord, 0, py );
908 GEOSGeometry* geosPt = GEOSGeom_createPoint_r( geosctxt, coord );
910 int type = GEOSGeomTypeId_r( geosctxt,
mGeos );
911 const GEOSGeometry* extRing = 0;
912 if ( type != GEOS_POLYGON )
919 extRing = GEOSGetExteriorRing_r( geosctxt,
mGeos );
921 GEOSCoordSequence *nearestCoord = GEOSNearestPoints_r( geosctxt, extRing, geosPt );
924 ( void )GEOSCoordSeq_getX_r( geosctxt, nearestCoord, 0, &nx );
925 ( void )GEOSCoordSeq_getY_r( geosctxt, nearestCoord, 0, &ny );
926 GEOSCoordSeq_destroy_r( geosctxt, nearestCoord );
927 GEOSGeom_destroy_r( geosctxt, geosPt );
936 catch ( GEOSException &e )
954 GEOSGeometry *centroidGeom = GEOSGetCentroid_r( geosctxt,
mGeos );
957 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, centroidGeom );
958 GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
959 GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
965 GEOSGeometry *pointGeom = GEOSPointOnSurface_r( geosctxt,
mGeos );
969 const GEOSCoordSequence *coordSeq = GEOSGeom_getCoordSeq_r( geosctxt, pointGeom );
970 GEOSCoordSeq_getX_r( geosctxt, coordSeq, 0, &px );
971 GEOSCoordSeq_getY_r( geosctxt, coordSeq, 0, &py );
972 GEOSGeom_destroy_r( geosctxt, pointGeom );
976 GEOSGeom_destroy_r( geosctxt, centroidGeom );
978 catch ( GEOSException &e )
994 while ( i <
nbPoints && ad[i] <= dl ) i++;
1004 di = sqrt( dx * dx + dy * dy );
1014 *px =
x[i] + dx * distr / di;
1015 *py =
y[i] + dy * distr / di;
1045 ( void )GEOSLength_r( geosctxt,
mGeos, &len );
1048 catch ( GEOSException &e )
const GEOSGeometry * geos() const
Returns the point set's GEOS geometry.
void getCentroid(double &px, double &py, bool forceInside=false) const
bool computeLineIntersection(double x1, double y1, double x2, double y2, double x3, double y3, double x4, double y4, double *x, double *y)
double dist_euc2d_sq(double x1, double y1, double x2, double y2)
struct pal::_cHullBox CHullBox
CHullBox * compute_chull_bbox()
QString tr(const char *sourceText, const char *disambiguation, int n)
bool qgsDoubleNear(double a, double b, double epsilon=4 *DBL_EPSILON)
void createGeosGeom() const
bool containsLabelCandidate(double x, double y, double width, double height, double alpha=0) const
Tests whether a possible label candidate will fit completely within the shape.
static void logMessage(const QString &message, const QString &tag=QString::null, MessageLevel level=WARNING)
add a message to the instance (and create it if necessary)
double dist_euc2d(double x1, double y1, double x2, double y2)
GEOSContextHandle_t geosContext()
Get GEOS context handle to be used in all GEOS library calls with reentrant API.
double cross_product(double x1, double y1, double x2, double y2, double x3, double y3)
bool containsPoint(double x, double y) const
Tests whether point set contains a specified point.
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)
double minDistanceToPoint(double px, double py, double *rx=0, double *ry=0) const
Returns the squared minimum distance between the point set geometry and the point (px...
const GEOSPreparedGeometry * preparedGeom() const
double length() const
Returns length of line geometry.
void getPointByDistance(double *d, double *ad, double dl, double *px, double *py)
Get a point a set distance along a line geometry.
static void splitPolygons(QLinkedList< PointSet * > &shapes_toProcess, QLinkedList< PointSet * > &shapes_final, double xrm, double yrm, const QgsFeatureId &uid)
Split a concave shape into several convex shapes.
void append(const T &value)
PointSet * extractShape(int nbPtSh, int imin, int imax, int fps, int fpe, double fptx, double fpty)