How to apply a method of interpolation (Voronoi polygons)?

Topics: SharpMap v0.9 / v1.x, SharpMap v2.0
Dec 11, 2007 at 7:54 PM
Is there any way to make it with SharpMap?
I can display a map from a point based shapefile with SharpMap, but I want to create a Contour layer and a GRID layer from a point based Layer.
Anyone can help me?
Dec 12, 2007 at 10:05 AM
These kinds of operations to not exist in SharpMap (and personally, I don't think they should). We've done some surfacing routines though, and when GeoAPI and NTS get fully refactored, we'll probably post these and some other full layer functions that we've done and make them available as a geo-processing support library. Probably won't happen until January, though.
Dec 12, 2007 at 12:07 PM
Thank Magnum4610, I will wait for these.
Feb 10, 2009 at 7:19 PM
Edited Feb 10, 2009 at 7:20 PM
Sharp map is a great software!
Well i am also intrested in this topic.
I also don't think that these kind of operations should be included within sharpmap, as they propably make it heavy.
If anyone has done some work i would be glad to hear from him/her!
I want to create either voronoi polygons or any another method to
create  regions from points.
Regards,
Agelos

Feb 15, 2009 at 1:35 PM
there is a voronoi clustering algorithm in C#  implemented (see http://bdittes.googlepages.com/ )
and some articles on Codeproject.com

Fortune's Voronoi algorithm implemented in C#
http://www.codeproject.com/KB/recipes/fortunevoronoi.aspx
Create temperature maps with 2D Voronoi diagrams
http://www.codeproject.com/KB/graphics/TemperatureVoronoi.aspx
 Visualization of the 2D Voronoi Diagram and the Delaunay Triangulation
http://www.codeproject.com/KB/graphics/VoronoiVisualization.aspx

Any ideas on how to use this code with SharpMap.Geometries.Point ?

I am having problems with comparing two points as it uses doubles as key for hash tables.

All the best,
Agelos

Feb 15, 2009 at 1:49 PM

I don’t have time to look further, but the Vector object in the project forces doubles to 10 digits of precision for the hash.   Can you overload the generator class to accept a list of SharpMap points, and cause it to convert the sharp map geometries to 10 digit representations?

From: agelospanagiotaki [mailto:notifications@codeplex.com]
Sent: Sunday, February 15, 2009 9:35 AM
To: tmacy@mapshots.com
Subject: Re: How to apply a method of interpolation (Voronoi polygons)? [SharpMap:18994]

From: agelospanagiotaki

there is a voronoi clustering algorithm in C# implemented (see http://bdittes.googlepages.com/ )
and some articles on Codeproject.com

Fortune's Voronoi algorithm implemented in C#
http://www.codeproject.com/KB/recipes/fortunevoronoi.aspx
Create temperature maps with 2D Voronoi diagrams
http://www.codeproject.com/KB/graphics/TemperatureVoronoi.aspx
Visualization of the 2D Voronoi Diagram and the Delaunay Triangulation
http://www.codeproject.com/KB/graphics/VoronoiVisualization.aspx

Any ideas on how to use this code with SharpMap.Geometries.Point ?

I am having problems with comparing two points as it uses doubles as key for hash tables.

All the best,
Agelos

Read the full discussion online.

To add a post to this discussion, reply to this email (SharpMap@discussions.codeplex.com)

To start a new discussion for this project, email SharpMap@discussions.codeplex.com

You are receiving this email because you subscribed to this discussion on CodePlex. You can unsubscribe on codePlex.com.

Please note: Images and attachments will be removed from emails. Any posts to this discussion will also be available online at codeplex.com

Coordinator
Feb 15, 2009 at 6:07 PM
Edited Feb 15, 2009 at 6:46 PM
Hi Agelo, I had a quick look, how about something like this.. note: it returns an IEnumerable<Geometry> but each geometry is actually a LineString not a set of closed LinearRings or polygons -  I haven't validated it visually but it builds and runs ok - you will need to add handling for points at infinity.. hth jd
 

using System.Collections.Generic;
using System.Linq;
using BenTools.Mathematics;
using SharpMap.Data;
using SharpMap.Geometries;

namespace SharpMapVoronoi
{
    public static class VoronoiExtensions
    {
        public static VoronoiGraph ComputeVoronoiGraph(IEnumerable<FeatureDataRow> points)
        {
            return ComputeVoronoiGraph(points.Select(o => o.Geometry as Point));
        }

        public static VoronoiGraph ComputeVoronoiGraph(IEnumerable<Point> points)
        {
            return Fortune.ComputeVoronoiGraph(points.Where(o => o != null).Select(o => new Vector(new[] { o.X, o.Y }){Tag = o}));
        }

        public static IEnumerable<Geometry> GeometriesFromVoronoiGraph(VoronoiGraph graph)
        {
            foreach (VoronoiEdge e in graph.Edges)
                yield return
                    new LineString(new List<double[]> { new[] { e.VVertexA[0], e.VVertexA[1] }, new[] { e.VVertexB[0], e.VVertexB[1] } });
        }

        public static IEnumerable<Geometry> VoronoiGeometriesFromFeatureDataTable(FeatureDataTable pointsTable)
        {
            return GeometriesFromVoronoiGraph(ComputeVoronoiGraph(pointsTable.AsEnumerable()));
        }

        public static IEnumerable<Geometry> VoronoiGeometriesFromPoints(IEnumerable<Point>  points)
        {
            return GeometriesFromVoronoiGraph(ComputeVoronoiGraph(points));
        }

        public static IEnumerable<FeatureDataRow> AsEnumerable(this FeatureDataTable dt)
        {
            foreach (FeatureDataRow r in dt.Rows)
                yield return r;
        }

    }
}


test 

using System.Diagnostics;
using Microsoft.VisualStudio.TestTools.UnitTesting;
using SharpMap.Data;
using SharpMap.Data.Providers;
using SharpMap.Geometries;

namespace SharpMapVoronoi.Test
{
    /// <summary>
    /// Summary description for TestSharpMapToVoronoi
    /// </summary>
    [TestClass]
    public class TestSharpMapToVoronoi
    {
        /// <summary>
        ///Gets or sets the test context which provides
        ///information about and functionality for the current test run.
        ///</summary>
        public TestContext TestContext { get; set; }

        #region Additional test attributes

        //
        // You can use the following additional attributes as you write your tests:
        //
        // Use ClassInitialize to run code before running the first test in the class
        // [ClassInitialize()]
        // public static void MyClassInitialize(TestContext testContext) { }
        //
        // Use ClassCleanup to run code after all tests in a class have run
        // [ClassCleanup()]
        // public static void MyClassCleanup() { }
        //
        // Use TestInitialize to run code before running each test
        // [TestInitialize()]
        // public void MyTestInitialize() { }
        //
        // Use TestCleanup to run code after each test has run
        // [TestCleanup()]
        // public void MyTestCleanup() { }
        //

        #endregion

        [TestMethod]
        public void TestMethod1()
        {
            var sf = new ShapeFile(@"xxxxxxxxx.shp");
            sf.Open();
            var fds = new FeatureDataSet();
            sf.ExecuteIntersectionQuery(sf.GetExtents(), fds);
            foreach (Geometry g in VoronoiExtensions.VoronoiGeometriesFromFeatureDataTable(fds.Tables[0]))
            {
                Debug.WriteLine(g.ToString());
            }
        }
    }
}

Feb 15, 2009 at 9:10 PM
Thank you very much for the help. 

I am unfortunately not using c# like you do, but anyway i found my way through. I got it to compile using .net 3 for the linq part.

I get an unhandled exception of type 'System.NotSupportedException' for my loaded point shape file.
dt.Rows seems to be the problem in this function.
 public static IEnumerable<FeatureDataRow> AsEnumerable(this FeatureDataTable dt)
        {
            if (dt.Rows != null)
            {
                foreach (FeatureDataRow r in dt.Rows)
                    yield return r;
            }
        }

any ideas?


Feb 15, 2009 at 9:38 PM
ok i changed it to
 foreach (FeatureDataRow r in dt)
and works ok.
THe problem is that the result is a hell!
I will investigate furter the issue and i will report here. I can upload also a sample if anyone is intrested.
Feb 15, 2009 at 9:49 PM
i uploaded a version here so you can check the differences
http://147.102.85.132/testsForVoronoi.rar

Coordinator
Feb 15, 2009 at 10:33 PM
Hi agelo, I progressed it a bit further, it is not solid but it will return polygons (though they are not guaranteed to be valid) and it should work with net 2 ;). It needs further work as it silently drops groups of edges which it cannot orientate hth jd 

using System;
using System.Collections.Generic;
using System.Collections.ObjectModel;
using BenTools.Mathematics;
using SharpMap.Data;
using SharpMap.Geometries;

namespace SharpMapVoronoi
{
    public static class VoronoiExtensions
    {
        public static VoronoiGraph ComputeVoronoiGraph(IEnumerable<FeatureDataRow> points)
        {
            return ComputeVoronoiGraph(L.Select(points, delegate(FeatureDataRow o) { return o.Geometry as Point; }));
        }

        public static VoronoiGraph ComputeVoronoiGraph(IEnumerable<Point> points)
        {
            return
                Fortune.ComputeVoronoiGraph(
                    L.Select(
                        L.Where(
                            points
                            , delegate(Point o) { return o != null; })
                        , delegate(Point o) { return new Vector(new[] { o.X, o.Y }) { Tag = o }; }));
        }

        public static IEnumerable<Geometry> GeometriesFromVoronoiGraph(VoronoiGraph graph, IEnumerable<Point> dataPoints)
        {
            IList<Point> points = L.ToList(dataPoints);

            Dictionary<Point, List<VoronoiEdge>> storage = new Dictionary<Point, List<VoronoiEdge>>();

            foreach (Point p in points)
            {
                foreach (VoronoiEdge e in graph.Edges)
                    if (e.LeftData.Tag == p || e.RightData.Tag == p)
                    {
                        if (!storage.ContainsKey(p))
                            storage.Add(p, new List<VoronoiEdge>());
                        storage[p].Add(e);
                    }
            }

            foreach (List<VoronoiEdge> edges in storage.Values)
            {
                Geometry g = EdgesToPolygon(edges);
                if (g != null)
                    yield return g;
            }
        }

        private static Geometry EdgesToPolygon(List<VoronoiEdge> edges)
        {
            Collection<Point> pts = new Collection<Point>();
            IEnumerable<VoronoiEdge> orientatedEdges = OrientateEdges(edges);
            if (orientatedEdges == null)
                return null;

            foreach (VoronoiEdge edge in orientatedEdges)
            {
                if (pts.Count == 0)
                    pts.Add(new Point(edge.VVertexA[0], edge.VVertexA[1]));
                pts.Add(new Point(edge.VVertexB[0], edge.VVertexB[1]));
            }
            return new Polygon(new LinearRing(pts));
        }
        //TODO:This method is incomplete and silently drops groups which it cannot orientate
        private static IEnumerable<VoronoiEdge> OrientateEdges(IEnumerable<VoronoiEdge> edges)
        {
            IList<VoronoiEdge> workingSet = new List<VoronoiEdge>(edges);
            IList<VoronoiEdge> result = new List<VoronoiEdge>();

            result.Add(workingSet[0]);
            workingSet.RemoveAt(0);

            while (workingSet.Count > 0)
            {
                int count = workingSet.Count;

                for (int i = workingSet.Count - 1; i > -1; i--)
                {
                    VoronoiEdge curr = workingSet[i];

                    bool found = false;
                    if (Equals(L.Last(result).VVertexB, curr.VVertexA))
                    {
                        found = true;
                        result.Add(workingSet[i]);
                    }
                    else if (Equals(L.Last(result).VVertexB, curr.VVertexB))
                    {
                        found = true;
                        result.Add(new VoronoiEdge { VVertexA = workingSet[i].VVertexB, VVertexB = workingSet[i].VVertexA });
                    }
                    else if (Equals(L.First(result).VVertexA, curr.VVertexB))
                    {
                        found = true;
                        result.Insert(0, workingSet[i]);
                    }
                    else if (Equals(L.First(result).VVertexA, curr.VVertexA))
                    {
                        found = true;
                        result.Insert(0,
                                      new VoronoiEdge { VVertexA = workingSet[i].VVertexB, VVertexB = workingSet[i].VVertexA });
                    }

                    if (found)
                        workingSet.RemoveAt(i);
                }

                if (workingSet.Count == count)//we didn't manage to orientate the list
                {
                    //throw new ApplicationException("Unable to orientate edges");//throw exception
                    return null;//silently drop
                }
            }
            return result;
        }

        public static bool Equals(Vector a, Vector b)
        {
            if (a.Dim == b.Dim)
            {
                for (int i = 0; i < a.Dim; i++)
                    if (a[i] != b[i])
                        return false;
                return true;
            }

            return false;
        }

        public static IEnumerable<Geometry> VoronoiGeometriesFromFeatureDataTable(FeatureDataTable pointsTable)
        {
            return
                VoronoiGeometriesFromPoints(
                    L.Select(L.EnumerateFeatureDataRows(pointsTable),
                             delegate(FeatureDataRow o) { return o.Geometry as Point; }));
        }

        public static IEnumerable<Geometry> VoronoiGeometriesFromPoints(IEnumerable<Point> points)
        {
            return GeometriesFromVoronoiGraph(ComputeVoronoiGraph(points), points);
        }
    }

    public static class L
    {
        public delegate U LFunc<T, U>(T o);

        public static T First<T>(IEnumerable<T> input)
        {
            foreach (T val in input)
            {
                return val;
            }
            return default(T);
        }

        public static T Last<T>(IEnumerable<T> input)
        {
            var lst = input as IList<T>;
            if (lst != null)
            {
                if (lst.Count > 0)
                    return lst[lst.Count - 1];
            }

            T curr = default(T);
            foreach (T c in input)
                curr = c;
            return curr;
        }

        public static IList<T> ToList<T>(IEnumerable<T> input)
        {
            IList<T> ls = new List<T>();
            foreach (T li in input)
                ls.Add(li);

            return ls;
        }

        public static IEnumerable<T> Where<T>(IEnumerable<T> input, LFunc<T, bool> predicate)
        {
            foreach (T val in input)
                if (predicate(val))
                    yield return val;
        }

        public static IEnumerable<U> Select<T, U>(IEnumerable<T> input, LFunc<T, U> conv)
        {
            foreach (T val in input)
                yield return conv(val);
        }

        public static IEnumerable<FeatureDataRow> EnumerateFeatureDataRows(FeatureDataTable dt)
        {
            foreach (FeatureDataRow r in dt.Rows)
                yield return r;
        }
    }
}

Feb 21, 2009 at 3:03 PM
I think i am geting better using sharpmap.i think i made progress.
 i found an Optimized implementation of Delaunay triangulation algorithm by Paul Bourke 
herehttp://local.wasp.uwa.edu.au/~pbourke/papers/triangulate/index.html 
and ofcourse tryed to transfer it to vb.net and apply the algorithm to sharpmap. 
I uploaded a new version ( http://147.102.85.132/testsForVoronoi.rar ) so everyone intersted can have a look.
As you can see for some reason produces less trinagles than expected.
Anyway my intention is to finaly produce the voronoi polygons.
Regards,
Agelos
Feb 23, 2009 at 9:36 PM

I took the previous version offline.I uploaded a new one with a "non correct" triangulation as it currectly converts doubles to singles(any clues on how to deal with this issue properly?).
i am planning to move on and create those voronoi polygons (even if the triangles are not so precise)
So i want to create these triangles as polygons. Then compute the centroid of each such triangle/polygon.
I don't have a clue on how to create polygons from these series of triangles.Again any help will be very much appreciated.

Regards,
Agelos

Feb 24, 2009 at 6:43 AM
Hi Agelos

SharpMap is alwayes in debt to it's original founder Morten Nielsen. He has done a smal implementation of Poul Bourkes work you can have a look at. While he has not used shapmap as the basis, it should be pretty straight forward to use the supplied code with sharpMap. (use trinangle vertices to generate SharpMap polygon).

Find it here: http://www.sharpgis.net/post/2006/03/09/Delaunay-Triangulation-in-NET-20.aspx



Best Peter
Mar 1, 2009 at 2:00 PM
Edited Mar 1, 2009 at 2:39 PM
Searching through the internet i found i found Emgu CV (openCV a cross platform .Net wrapper to the Intel OpenCV image-processing library).
I created a demo with SharpMap which seems to correctly calculate both Delauny and Voronoi. 
Ofcourse unfortunately i convert Sharpmap.geometries.Point to System.Drawing.PointF. 
Is there any method to compare the differences of the produced voronoi and  with The result given?  I incude in the demo a layer of correctly calculated voronois.
I also want to chop the produced polygons against another layer(The country layer). Any ideas on how to do this?
Please download from here http://energy.chemeng.ntua.gr/openCV_Sharpmap_Voronoi_Delauny.rar and advice accordingly.

Best Regards,
Agelos
Mar 2, 2009 at 5:31 PM
Edited Mar 2, 2009 at 5:32 PM
I have managed to chop the voronoi polygons based on the country layer. And also create a buffer of it using NTS. 
I uploaded a new version for you to see:
Coordinator
Mar 2, 2009 at 9:09 PM
Looking good Agelos!
Mar 19, 2009 at 10:42 AM
Hi all, first of all, thanks for developing this great and very useful tool!
I have an exception when click on calculate button (on the choped version) about the init of Emgu.CV.CvInvoke.

Can you help me please??

Thanks at all in advance!!
Dado

Mar 19, 2009 at 12:44 PM
Can you please post the exception here?
Mar 19, 2009 at 1:23 PM
Hi,
the exception occour on the following row

Dim subdivision As PlanarSubdivision = New PlanarSubdivision(verticesPF)

of the sub CalculateBothDelaunyVoronoi() in the frmMain.

it is:

An unhandled exception of type 'System.TypeInitializationException' occurred in Emgu.CV.dll
Additional information: L'inizializzatore di tipo di 'Emgu.CV.CvInvoke' ha generato un'eccezione.

(sorry, it's in italian..)

Thanks,
Dado
Mar 19, 2009 at 4:25 PM
I can imagine what the problem is, even though .... i should learn italian! :)
As this program requires Emgu.CV  and Open CV just copy all the dlls provided by Emgu.CV ( http://emgu.com/wiki/index.php) in the folder of /debug? or release.  
hth otherwise please contact me again.
Mar 23, 2009 at 3:03 PM
Hi,
I've solved my problem: I have copied all DLLs (as you suggested) but it was not enough.
I've installed also C++ Redistributable package to solve completely my problem.

Now, all works, thanksa lot!!!

Bye,
Dado