← Back to branch summary

~sharpie/geos/3.3.2

« back to all changes in this revision

Viewing changes to source/algorithm/RobustLineIntersector.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Fabio Tranchitella
  • Date: 2006-11-06 21:35:52 UTC
  • mfrom: (3.1.3 feisty)
  • Revision ID: james.westby@ubuntu.com-20061106213552-m03o92ggj1na737b
Tags: 2.2.3-3
debian/control: move doxygen from build-depends-indep to build-depends.

Show diffs side-by-side

added added

removed removed

Lines of Context:
source/algorithm/RobustLineIntersector.cpp
1
1
/**********************************************************************
2
 
 * $Id: RobustLineIntersector.cpp,v 1.29 2004/12/08 13:54:43 strk Exp $
 
2
 * $Id: RobustLineIntersector.cpp,v 1.29.2.1 2005/05/23 16:39:37 strk Exp $
3
3
 *
4
4
 * GEOS - Geometry Engine Open Source
5
5
 * http://geos.refractions.net
6
6
 *
7
7
 * Copyright (C) 2001-2002 Vivid Solutions Inc.
 
8
 * Copyright (C) 2005 Refractions Research Inc.
8
9
 *
9
10
 * This is free software; you can redistribute and/or modify it under
10
11
 * the terms of the GNU Lesser General Public Licence as published
22
23
#define COMPUTE_Z 1
23
24
#endif // COMPUTE_Z
24
25
 
 
26
namespace { // module-statics
 
27
 
 
28
using namespace geos;
 
29
 
 
30
void
 
31
normalizeToEnvCentre(Coordinate &n00, Coordinate &n01,
 
32
                Coordinate &n10, Coordinate &n11, Coordinate &normPt) 
 
33
{
 
34
        double minX0 = n00.x < n01.x ? n00.x : n01.x;
 
35
        double minY0 = n00.y < n01.y ? n00.y : n01.y;
 
36
        double maxX0 = n00.x > n01.x ? n00.x : n01.x;
 
37
        double maxY0 = n00.y > n01.y ? n00.y : n01.y;
 
38
 
 
39
        double minX1 = n10.x < n11.x ? n10.x : n11.x;
 
40
        double minY1 = n10.y < n11.y ? n10.y : n11.y;
 
41
        double maxX1 = n10.x > n11.x ? n10.x : n11.x;
 
42
        double maxY1 = n10.y > n11.y ? n10.y : n11.y;
 
43
 
 
44
        double intMinX = minX0 > minX1 ? minX0 : minX1;
 
45
        double intMaxX = maxX0 < maxX1 ? maxX0 : maxX1;
 
46
        double intMinY = minY0 > minY1 ? minY0 : minY1;
 
47
        double intMaxY = maxY0 < maxY1 ? maxY0 : maxY1;
 
48
 
 
49
        double intMidX = (intMinX + intMaxX) / 2.0;
 
50
        double intMidY = (intMinY + intMaxY) / 2.0;
 
51
 
 
52
        normPt.x = intMidX;
 
53
        normPt.y = intMidY;
 
54
 
 
55
        n00.x -= normPt.x;    n00.y -= normPt.y;
 
56
        n01.x -= normPt.x;    n01.y -= normPt.y;
 
57
        n10.x -= normPt.x;    n10.y -= normPt.y;
 
58
        n11.x -= normPt.x;    n11.y -= normPt.y;
 
59
 
 
60
#if COMPUTE_Z
 
61
        double minZ0 = n00.z < n01.z ? n00.z : n01.z;
 
62
        double minZ1 = n10.z < n11.z ? n10.z : n11.z;
 
63
        double maxZ0 = n00.z > n01.z ? n00.z : n01.z;
 
64
        double maxZ1 = n10.z > n11.z ? n10.z : n11.z;
 
65
        double intMinZ = minZ0 > minZ1 ? minZ0 : minZ1;
 
66
        double intMaxZ = maxZ0 < maxZ1 ? maxZ0 : maxZ1;
 
67
        double intMidZ = (intMinZ + intMaxZ) / 2.0;
 
68
        normPt.z = intMidZ;
 
69
        n00.z -= normPt.z;
 
70
        n01.z -= normPt.z;
 
71
        n10.z -= normPt.z;
 
72
        n11.z -= normPt.z;
 
73
#endif
 
74
}
 
75
 
 
76
} // module-statics
 
77
 
25
78
namespace geos {
26
79
 
27
80
RobustLineIntersector::RobustLineIntersector(){}
385
438
Coordinate*
386
439
RobustLineIntersector::intersection(const Coordinate& p1,const Coordinate& p2,const Coordinate& q1,const Coordinate& q2) const
387
440
{
388
 
        Coordinate *n1=new Coordinate(p1);
389
 
        Coordinate *n2=new Coordinate(p2);
390
 
        Coordinate *n3=new Coordinate(q1);
391
 
        Coordinate *n4=new Coordinate(q2);
392
 
        Coordinate *normPt=new Coordinate();
393
 
        normalize(n1, n2, n3, n4, normPt);
394
 
        Coordinate intPt;
 
441
        Coordinate n1=p1;
 
442
        Coordinate n2=p2;
 
443
        Coordinate n3=q1;
 
444
        Coordinate n4=q2;
 
445
        Coordinate normPt;
 
446
 
 
447
        //normalize(n1, n2, n3, n4, normPt);
 
448
        normalizeToEnvCentre(n1, n2, n3, n4, normPt);
 
449
 
 
450
        Coordinate *intPt=NULL;
 
451
 
395
452
        try {
396
 
                Coordinate *h=HCoordinate::intersection(*n1,*n2,*n3,*n4);
397
 
                intPt.setCoordinate(*h);
398
 
#if COMPUTE_Z
399
 
                double ztot = 0;
400
 
                double zvals = 0;
401
 
                double zp = interpolateZ(intPt, p1, p2);
402
 
                double zq = interpolateZ(intPt, q1, q2);
403
 
                if ( !ISNAN(zp)) { ztot += zp; zvals++; }
404
 
                if ( !ISNAN(zq)) { ztot += zq; zvals++; }
405
 
                if ( zvals ) intPt.z = ztot/zvals;
406
 
#endif // COMPUTE_Z
407
 
                delete h;
 
453
                intPt=HCoordinate::intersection(n1,n2,n3,n4);
408
454
        } catch (NotRepresentableException *e) {
 
455
                delete intPt;
409
456
                Assert::shouldNeverReachHere("Coordinate for intersection is not calculable"+e->toString());
410
 
    }
411
 
        intPt.x+=normPt->x;
412
 
        intPt.y+=normPt->y;
413
 
        if (precisionModel!=NULL) {
414
 
                precisionModel->makePrecise(&intPt);
415
457
        }
416
 
        delete n1;
417
 
        delete n2;
418
 
        delete n3;
419
 
        delete n4;
420
 
        delete normPt;
421
 
    /**
422
 
     * MD - after fairly extensive testing
423
 
     * it appears that the computed intPt always lies in the segment envelopes
424
 
     */
425
 
    //if (! isInSegmentEnvelopes(intPt))
426
 
    //    System.out.println("outside segment envelopes: " + intPt);
427
 
 
428
 
        return new Coordinate(intPt);
 
458
 
 
459
        intPt->x += normPt.x;
 
460
        intPt->y += normPt.y;
 
461
 
 
462
/**
 
463
 *
 
464
 * MD - May 4 2005 - This is still a problem.  Here is a failure case:
 
465
 *
 
466
 * LINESTRING (2089426.5233462777 1180182.3877339689,
 
467
 *             2085646.6891757075 1195618.7333999649)
 
468
 * LINESTRING (1889281.8148903656 1997547.0560044837,
 
469
 *             2259977.3672235999 483675.17050843034)
 
470
 * int point = (2097408.2633752143,1144595.8008114607)
 
471
 */
 
472
 
 
473
        if (precisionModel!=NULL) precisionModel->makePrecise(intPt);
 
474
 
 
475
 
 
476
#if COMPUTE_Z
 
477
        double ztot = 0;
 
478
        double zvals = 0;
 
479
        double zp = interpolateZ(*intPt, p1, p2);
 
480
        double zq = interpolateZ(*intPt, q1, q2);
 
481
        if ( !ISNAN(zp)) { ztot += zp; zvals++; }
 
482
        if ( !ISNAN(zq)) { ztot += zq; zvals++; }
 
483
        if ( zvals ) intPt->z = ztot/zvals;
 
484
#endif // COMPUTE_Z
 
485
 
 
486
        return intPt;
429
487
}
430
488
 
431
489
 
494
552
 
495
553
/**********************************************************************
496
554
 * $Log: RobustLineIntersector.cpp,v $
 
555
 * Revision 1.29.2.1  2005/05/23 16:39:37  strk
 
556
 * Ported JTS robustness patch for RobustLineIntersector
 
557
 *
497
558
 * Revision 1.29  2004/12/08 13:54:43  strk
498
559
 * gcc warnings checked and fixed, general cleanups.
499
560
 *