STIntersection performance

[FredL]

Hi,

I've been getting some slow responses from STIntersection() using NTS,
and have tried performing the same intersection with Postgis and
Microsoft.SqlServer.Types.dll to see how they compare. The differences
are remarkable:

Postgis: 300 - 400ms
SqlServer.Types.dll: 50 - 100ms
Nts: 2 minutes, 21 seconds!

I'm using version 1.7.3 on an XP VM. (I also tried running the NTS
intersection test using mono on a mac with similar results). I've
tried to add the WKT strings for the two shapes to this message, but
the post was rejected, probably because one of the strings was too
large - will see if I can attach.

Any suggestions would be most welcome!

Many thanks.

[FredL]

WKT strings attached.

[Diego Guidi]

I twy to take a look as fast as I can, checkin how fast is also JTS
suite doing same task

Diego Guidi

On Sun, Apr 24, 2011 at 1:15 PM, FredL <fred.l...@gmail.com> wrote:
> WKT strings attached.
>
> --
> You received this message because you are subscribed to the Google Groups
> "NetTopologySuite" group.
> To post to this group, send email to nettopol...@googlegroups.com.
> To unsubscribe from this group, send email to
> nettopologysui...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/nettopologysuite?hl=en.
>

[FredL]

Thanks Diego.

On Apr 26, 7:25 am, Diego Guidi <diegogu...@gmail.com> wrote:
> I twy to take a look as fast as I can, checkin how fast is also JTS
> suite doing same task
>
> Diego Guidi
>

[Diego Guidi]

the same test runs in about ten seconds in JTS 1.8... now it's time to find how NTS loses all the time

[FredL]

thanks again. Let me know if I can help at this end...

[Felix Obermaier]

The same test takes about 1.111ms using NTS v2.11, so I think it is related to changes in the algorithm.
Cheers FObermaier

> -----Ursprüngliche Nachricht-----
> Von: nettopol...@googlegroups.com
> [mailto:nettopol...@googlegroups.com] Im Auftrag von FredL
> Gesendet: Dienstag, 26. April 2011 15:14
> An: NetTopologySuite
> Betreff: Re: R: Re: STIntersection performance

[FredL]

Hi Felix,

Could you let me know what the svn branch is for v2.11? I'll try a
build at this end and test performance..

[Diego Guidi]

> The same test takes about 1.111ms using NTS v2.11, so I think it is related to changes in the algorithm.

I didn't checked nothing at all, but I suspect that can be related to
collections used by NTS, i.e. Set/HashSet, PowerCollections, and so
on... anything related to this one?
http://code.google.com/p/nettopologysuite/issues/detail?id=43&can=1
this issue is fixed, at least in 2.11

[Diego Guidi]

OverlayOp.CopyPoints for the multipolygon geometry is executed in a fraction of a second by JTS, and in 10 seconds or more by NTS (trunk)

I think that here is the main problem, now the hard part: discover why this happens :)

[John Diss]

Hi Diego, perhaps try with r538 or earlier to make sure I didn’t break anything aliasing ArrayList = List<Object> etc when doing the Silverlight stuff...

--

[Diego Guidi]

sorry, I'm wrong. Maybe a temporary result due to some problem in my pc... I check better and see

[John Diss]

 

Thinking about it the alias only applies within Silverlight – in .Net it stays as ArrayList...

[diego...@gmail.com]

If I disable NodingValidator (OverlayOp.NodingValidatorDisabled = true;) the operation is executed immediately.
Of course this is not a good thing to do, and JTS validates noding, so I think that is the nodingvalidator that actually is slow.

Il giorno , Diego Guidi <diego...@gmail.com> ha scritto:

> sorry, I'm wrong. Maybe a temporary result due to some problem in my pc... I check better and see
>
>
>
>

[Diego Guidi]

> Thinking about it the alias only applies within Silverlight – in .Net it
> stays as ArrayList...

exactly

[John Diss]

I didn't break it :)
::dances::

-----Original Message-----
From: nettopol...@googlegroups.com [mailto:nettopol...@googlegroups.com] On Behalf Of Diego Guidi
Sent: 26 April 2011 16:42
To: nettopol...@googlegroups.com
Subject: Re: R: Re: STIntersection performance

> Thinking about it the alias only applies within Silverlight - in .Net it
> stays as ArrayList...
exactly

--

[Diego Guidi]

> I didn't break it :)
> ::dances::

NodingValidator.CheckInteriorIntersections looks the only responsible
of this poor performances, not you :)

[Diego Guidi]

CheckEndPtVertexIntersections 00:00:00.1641412
CheckInteriorIntersections 00:02:11.3173618
CheckCollapses 00:00:00.0017678

[FObermaier]

NodingValidator has been superceeded by FastNodingValidator
cheers Felix

On Apr 26, 5:48 pm, Diego Guidi <diegogu...@gmail.com> wrote:
> CheckEndPtVertexIntersections 00:00:00.1641412
> *CheckInteriorIntersections 00:02:11.3173618*
> CheckCollapses 00:00:00.0017678

[Diego Guidi]

> NodingValidator has been superceeded by FastNodingValidator

this isn't true for the trunk

[Diego Guidi]

FastNodingValidator is used by 2.x branch, and looks that contains "only" check for interior intersections, actually the same heck that is so slow in trunk. So maybe here's why trunk and 2.x branches are so slow.

[Diego Guidi]

I've added a test that validates a single segmentstring (the slower of all the checks, alone it executes in 5 seconds).

In NTS the check is executed in 5 seconds (as expected), in JTS the execution is immediate.

I hope that analyzing the single test I can find how is the problem

[Diego Guidi]

why god simply using a one-dimensional array can improve performances to 2x or more?

the test on segmentstring now executes in 2 seconds instead of 5, and all I've changed is that now in LineIntersector

protected ICoordinate[] inputLines = new ICoordinate[4];

and not 

protected ICoordinate[,] inputLines = new ICoordinate[2,2];

[FredL]

well found. It seems v.strange that a multi-dimensional array would be
so much slower...

[Diego Guidi]

> well found. It seems v.strange that a multi-dimensional array would be
> so much slower...

not the only reason to slow performances... actually I don't
understand why and how NTS is slow

[Diego Guidi]

> well found. It seems v.strange that a multi-dimensional array would be
> so much slower...

the test SlowIntersectionTest.TestAsDefault runs now:
in 26 seconds with NTS-Release
and in 10 seconds in JTS 1.8.0

this is the best at the moment :(

[Felix Obermaier]

Diego, I started backporting FastNodingValidator to trunk, maybe there is some more improvement there
Felix

> -----Ursprüngliche Nachricht-----
> Von: nettopol...@googlegroups.com
> [mailto:nettopol...@googlegroups.com] Im Auftrag von Diego Guidi
> Gesendet: Dienstag, 3. Mai 2011 09:07
> An: nettopol...@googlegroups.com
> Betreff: Re: R: Re: R: Re: R: Re: STIntersection performance

[Diego Guidi]

> Diego, I started backporting FastNodingValidator to trunk,
> maybe there is some more improvement there

thanks, anyway what i'm experiencing it's strange:
Look at the code:
LineIntersector.ComputeIntersection(ICoordinate p1, ICoordinate p2,
ICoordinate p3, ICoordinate p4)
if i use this code for the method:
inputLines[0] = p1;
inputLines[1] = p2;
inputLines[2] = p3;
inputLines[3] = p4;
//result = ComputeIntersect(p1, p2, p3, p4);

slows the test SlowIntersectionTest.TestSlowerSegmentString by 8 times
(from 100 to 800 milliseconds, in Release). but inputLines isn't used,
is simply filled!

[Diego Guidi]

> slows the test SlowIntersectionTest.TestSlowerSegmentString by 8 times
> (from 100 to 800 milliseconds, in Release). but inputLines isn't used,
> is simply filled!

times are comparable in other cases, i.e
//inputLines[0] = p1;
//inputLines[1] = p2;
//inputLines[2] = p3;
//inputLines[3] = p4;

result = ComputeIntersect(p1, p2, p3, p4);

executes in the same time in java and .net

[FObermaier]

how about
inputLines = new[] {p1, p2, p3, p4};

DotTrace is a really good profiler, that even works in NUnit tests.

[diego...@gmail.com]

> how about
> inputLines = new[] {p1, p2, p3, p4};

same performances, I can try better with a dottrace profiler trial or with another profiler

[Felix Obermaier]

Diego,

 

here some first results:

SlowIntersectionTest.TestSlowerSegmentString : Passed

NodingValidator.CheckValid => ElapsedMilliseconds: 1741
FastNodingValidator.CheckValid => ElapsedMilliseconds: 62

 

 

--

[Diego Guidi]

> NodingValidator.CheckValid => ElapsedMilliseconds: 1741
> FastNodingValidator.CheckValid => ElapsedMilliseconds: 62

wow :-0
Feel free to commit in the trunk, I check as soon as possible

[Diego Guidi]

thanks to work of Felix Obermaier (see: http://code.google.com/p/nettopologysuite/source/detail?r=566)

now FastNodingValidator is used instead of standard NodingValidator.

Now intersections are almost immediate. a really great improvement, thanks a lot Felix.

[FredL]

Almost a year on I thought I should say thanks for fixing this - will re-try implementing NTS in my project.

[Jeremy]

Is this on 1.12.1?  I'm seeing pretty awful times for Intersection.

[FObermaier]

It is in 1.12.1.
If you have to perform a one to many (1:n) check, use functionality of PreparedGeometry.

[Jeremy]

I don't see how to use PreparedGeometry to calculate ST_Intersection.  It only contains a call for Intersects.  I'm trying to do something like:

            // Create a grid of polygons over the input geometry's envelope.
            List<IGeometry> gridCells = new List<IGeometry>(numCols * numRows);

            for (int y = 0; y < numRows; y++)
            {
                for (int x = 0; x < numCols; x++ )
                {
                    ILinearRing ring = new LinearRing(new Coordinate[]
                    {
                        new Coordinate(xOffset + (x * cellWidth), yOffset + (y * cellHeight)),
                        new Coordinate(xOffset + ((x + 1) * cellWidth), yOffset + (y * cellHeight)),
                        new Coordinate(xOffset + ((x + 1) * cellWidth), yOffset + ((y + 1) * cellHeight)),
                        new Coordinate(xOffset + (x * cellWidth), yOffset + ((y + 1) * cellHeight)),
                        new Coordinate(xOffset + (x * cellWidth), yOffset + (y * cellHeight))
                    });

                    gridCells.Add(new Polygon(ring));
                }
            }

            IGeometry collection = new GeometryCollection(gridCells.ToArray());

            // Compute the intersection of the grid (removes the parts that are outside of the actual polygon)
            IGeometry intersection = collection.Intersection(shapeTransformed);

[Felix Obermaier]

You are right, PreparedGeometries are for predicate computation only.

 

To solve a problem like yours, you might want to

-          Create a spatial index (quadtree) on your grid geometries.

-          Query that for possible intersecting grid cells

-          Perform a intersects operation on only those.

-          Union the resulting geometries (if that is necessary at all)

 

[Jeremy]

Unfortunately I use the original geometry's envelope to calculate the offsets/widths/heights in the loop above so the quadtree doesn't cull anything out:

public static IGeometry GridSpan(Context context, IGeometry shape, double metersX, double metersY)
        {
            Stopwatch sw = new Stopwatch();
            sw.Start();
           
            uint zone;
            bool isSouthernHemisphere;
            Wgs84ToUtmParameters(shape.Coordinate.Y, shape.Coordinate.X, out zone, out isSouthernHemisphere);

            GeographicCoordinateSystem from = GeographicCoordinateSystem.WGS84;
            ProjectedCoordinateSystem to = ProjectedCoordinateSystem.WGS84_UTM((int)zone, !isSouthernHemisphere);

            IGeometry shapeTransformed = TransformGeometry(new CoordinateTransformationFactory().CreateFromCoordinateSystems(from, to), shape);

            IGeometry grid = shapeTransformed.Envelope;
            Coordinate[] gridCoordinates = grid.Coordinates; // Everytime the Coordinates property is called it generates a big array so cache it.

            double gridWidth = gridCoordinates[2].X - gridCoordinates[0].X;
            double gridHeight = gridCoordinates[2].Y - gridCoordinates[0].Y; // 0 was 1 on the web page?

            // Lower left
            double xOffset = gridCoordinates[0].X;
            double yOffset = gridCoordinates[0].Y;

            double cellWidth = metersX;
            double cellHeight = metersY;
            int numCols = (int)Math.Ceiling(gridWidth / cellWidth);
            int numRows = (int)Math.Ceiling(gridHeight / cellHeight);

            // Create the grid
            Quadtree<IPolygon> quadTree = new Quadtree<IPolygon>();

            for (int y = 0; y < numRows; y++)
            {
                for (int x = 0; x < numCols; x++ )
                {
                    ILinearRing ring = new LinearRing(new Coordinate[]
                    {
                        new Coordinate(xOffset + (x * cellWidth), yOffset + (y * cellHeight)),
                        new Coordinate(xOffset + ((x + 1) * cellWidth), yOffset + (y * cellHeight)),
                        new Coordinate(xOffset + ((x + 1) * cellWidth), yOffset + ((y + 1) * cellHeight)),
                        new Coordinate(xOffset + (x * cellWidth), yOffset + ((y + 1) * cellHeight)),
                        new Coordinate(xOffset + (x * cellWidth), yOffset + (y * cellHeight))
                    });

                    quadTree.Insert(new Envelope(xOffset + (x * cellWidth), xOffset + ((x + 1) * cellWidth),
                                                 yOffset + (y * cellHeight), yOffset + ((y + 1) * cellHeight)),
                                    new Polygon(ring));
                }
            }

            IList<IPolygon> polys = quadTree.Query(shapeTransformed.EnvelopeInternal);
            IPolygon[] array = new IPolygon[polys.Count];
            polys.CopyTo(array, 0);

            NetTopologySuite.Operation.Overlay.OverlayOp.NodingValidatorDisabled = true;

            IGeometry collection = new GeometryCollection(array);
            IGeometry intersection = collection.Intersection(shapeTransformed);
            IGeometry transform = TransformGeometry(new CoordinateTransformationFactory().CreateFromCoordinateSystems(to, from), intersection);

            sw.Stop();
            Console.WriteLine("Grid Total (ms) " + sw.ElapsedMilliseconds);

            return transform;

[FObermaier]

You may have better performance using this:

public static IGeometry GridSpan(Context context, IGeometry shape, double metersX, double metersY)
        {

            var sw = new Stopwatch();
            var swIntersection = new StopWatch();

            sw.Start();

            uint zone;
            bool isSouthernHemisphere;
            Wgs84ToUtmParameters(shape.

Coordinate.Y, shape.Coordinate.X, out zone, out isSouthernHemisphere);

            GeographicCoordinateSystem from = GeographicCoordinateSystem.WGS84;
            ProjectedCoordinateSystem to = ProjectedCoordinateSystem.WGS84_UTM((int)zone, !isSouthernHemisphere);

            IGeometry shapeTransformed = TransformGeometry(new CoordinateTransformationFactory().CreateFromCoordinateSystems(from, to), shape);

            var gridEnv = shapeTransformed.EnvelopeInternal;

            double gridWidth = gridEnv.Width;
            double gridHeight = gridEnv.Height; // 0 was 1 on the web page?

            // Lower left
            double xOffset = gridEnv.MinX;
            double yOffset = gridEnv.MinY;

            double cellWidth = metersX;
            double cellHeight = metersY;
            int numCols = (int)Math.Ceiling(gridWidth / cellWidth);
            int numRows = (int)Math.Ceiling(gridHeight / cellHeight);

            var gridList = new List<IGeometry>(numCols * numRows)
            var factory = shapeTransformed.Factory;

            // Create the grid

            for (int y = 0; y < numRows; y++)
            {
                for (int x = 0; x < numCols; x++ )
                {

                    var ring = factory.CreateLinearRing(new []

                    {
                        new Coordinate(xOffset + (x * cellWidth), yOffset + (y * cellHeight)),
                        new Coordinate(xOffset + ((x + 1) * cellWidth), yOffset + (y * cellHeight)),
                        new Coordinate(xOffset + ((x + 1) * cellWidth), yOffset + ((y + 1) * cellHeight)),
                        new Coordinate(xOffset + (x * cellWidth), yOffset + ((y + 1) * cellHeight)),
                        new Coordinate(xOffset + (x * cellWidth), yOffset + (y * cellHeight))
                    });

                    var poly = factory.CreatePolygon(ring, null);
                    swIntersection.Start();
                    var intersection = transformedShape.Intersection(poly);
                    swIntersection.Stop();
                    if (intersection != null) gridList.Add(intersection);

                }
            }

            var collection = factory.CreateGeometryCollection(gridList.ToArray());
            IGeometry transform = TransformGeometry(new CoordinateTransformationFactory().CreateFromCoordinateSystems(to, from), collection);

            sw.Stop();
            Console.WriteLine("Grid Total (ms) " + sw.ElapsedMilliseconds);

            Console.WriteLine("Intersection Total (ms) " + swIntersection.ElapsedMilliseconds);

            return transform;

}

Hth FObermaier

[FObermaier]

Depending on the complexity of your input geometry you may get better performance by performing a "contains" check using prepared geometry

var pg = new NetTopologySuite.Geometries.Prepared.PreparedGeometryFactory().Prepare(transformedShape);

for ...
   for ...
 

[FObermaier]

   {
       ...
       if (pg.Contains(poly))
           gridList.Add(poly);
       else if (pg.Intersects(poly))
       {
           var int = transformedShape.Intersection(poly);
           if (int != null)
               gridList.Add(int);
       }
   }

 

}

             

[Jeremy]

That did the trick!  Thank you so much.  For dividing a 5 acre plot into 1x1 meter cells the time dropped from 40 seconds down to 8 seconds.  One question: you switched some of my code to using factories... is that prefered over explicit instantiation?

Thanks,
Jeremy

[Felix Obermaier]

> That did the trick!  Thank you so much.  For dividing a 5 acre plot into 1x1 meter cells the time dropped from 40 seconds down to 8 seconds. 

You're welcome.

> One question: you switched some of my code to using factories... is that prefered over explicit instantiation?

Yes