[gdal-dev] Ogr: Is an OGR::Geometry::InteriorPoint() method ?

[Benjamin]

Hi,



Even Rouault wrote
> 
> To my knowledge there is no such method. If you don't have particular
> constraints on the point except that it must be at the interior of the
> polygon,
> you could take a random vertex of the boundary, compute the distance D to
> the
> next point, and then compute various points located on a circle centered
> on that
> vertex and of radius D/2 until the IsPointOnSurface() returns TRUE. Not
> very
> efficient, there are likely better methods...
> 

I don't want a totally random point.



Ari Jolma-2 wrote
> 
> On 08/16/2012 03:09 PM, Even Rouault wrote:
>> Selon Benjamin <benjamin.lux@>:
>>
>> > For example, with JTS (Java Topologie Suite) we can obtain an interior
>> point (0.5 2.5).
>> Do you know the name of the JTS method ? It might also exist in GEOS, the
>> C++
>> port of JTS that GDAL relies on to do geometry related computations.
> 
> http://geos.osgeo.org/doxygen/classgeos_1_1algorithm_1_1InteriorPointArea.html 
> describes geos::algorithm::InteriorPointArea Class, which can be used to 
> obtain an interior point of an area.
> 
> Ari
> 

I had a look to JTS sources and write again the InteriorPointArea in C#.
If someone want it I take it at the end of this email (and for others .NET
users who google it).

Thank to Even and Ari for theirs answers.


<code>
// A copie from JTS InteriorPointArea Class (Java) for C#
// See also : http://www.vividsolutions.com/jts/jtshome.htm

using OSGeo.OGR;
/// <summary>
/// Computes a point in the interior of an area geometry.
/// 
/// Algorithm :
/// 1/ Find the intersections between the geometry 
/// and the horizontal bisector of the area's envelope
/// 2/Pick the midpoint of the largest intersection 
/// (the intersections will be lines and points)
/// </summary>
/// <remarks>
/// Note: If a fixed precision model is used,
/// in some cases this method may return a point
/// which does not lie in the interior.
/// </remarks>
public class InteriorPointArea
{

    private static double Avg(double a, double b)
    {
        return (a + b) / 2.0;
    }

    //private Geometry _factory;
    private double[] _interiorPoint = null;
    private double _maxWidth = 0.0;

    public InteriorPointArea(Geometry g)
    {
        //_factory = g;
        Add(g);
    }

    public double[] GetInteriorPoint()
    {
        return _interiorPoint;
    }

    /// <summary>
    /// Tests the interior vertices (if any)
    /// defined by a linear Geometry for the best inside point.
    /// If a Geometry is not of dimension 1 it is not tested.
    /// </summary>
    /// geom the geometry to add.
    private void Add(Geometry geom)
    {
        if (geom.GetGeometryType() == wkbGeometryType.wkbPolygon)
        {
            AddPolygon(geom);
        }
        else if (geom.GetGeometryType() ==
wkbGeometryType.wkbGeometryCollection
            || geom.GetGeometryType() == wkbGeometryType.wkbMultiPolygon)
        {
            for (int i = 0; i < geom.GetGeometryCount(); i++)
            {
                Add(geom.GetGeometryRef(i));
            }
        }
    }



    /// <summary>
    ///  Finds a reasonable point at which to label a Geometry.
    /// @param geometry the geometry to analyze
    /// </summary>
    /// the geometry to analyze.
    /// <returns>the midpoint of the largest intersection between the
geometry and
    ///  a line halfway down its envelope</returns>
    public void AddPolygon(Geometry geometry)
    {
        Geometry bisector = HorizontalBisector(geometry);

        Geometry intersections = bisector.Intersection(geometry);
        Geometry widestIntersection = WidestGeometry(intersections);

        double width = widestIntersection.Length();
        if (_interiorPoint == null || width > _maxWidth)
        {
            _interiorPoint = CentreOfEnveloppe(widestIntersection);
            _maxWidth = width;
        }
    }



    /// <summary>
    /// return the widests geometry.
    /// </summary>
    /// The geometry.
    /// <returns>
    /// if geometry is a collection, the widest sub-geometry; otherwise,
    /// the geometry itself
    /// </returns>
    protected Geometry WidestGeometry(Geometry geometry)
    {
        if (geometry.GetGeometryType() != wkbGeometryType.wkbMultiLineString
&&
            geometry.GetGeometryType() !=
wkbGeometryType.wkbGeometryCollection)
        {
            return geometry;
        }

        if (geometry.IsEmpty())
        {
            return geometry;
        }

        Geometry widestGeometry = geometry.GetGeometryRef(0);

        for (int i = 1; i < geometry.GetGeometryCount(); i++)
        { //Start at 1
            if (geometry.GetGeometryRef(i).Length() >
widestGeometry.Length())
            {
                widestGeometry = geometry.GetGeometryRef(i);
            }
        }
        return widestGeometry;
    }

    protected Geometry HorizontalBisector(Geometry geometry)
    {

        Envelope envelope = new Envelope();
        geometry.GetEnvelope(envelope);

        // Assert: for areas, minx <> maxx
        double avgY = Avg(envelope.MinY, envelope.MaxY);
        Geometry result = new Geometry(wkbGeometryType.wkbLineString);
        result.AddPoint(envelope.MinX, avgY, 0);
        result.AddPoint(envelope.MaxX, avgY, 0);
        return result;
    }

    /// <summary>
    /// Returns the centre point of the envelope.
    /// </summary>
    /// envelope the envelope to analyze.
    /// <returns>Returns the centre point of the envelope.</returns>
    public double[] Centre(Envelope envelope)
    {
        return new double[] { Avg(envelope.MinX, envelope.MaxX),
Avg(envelope.MinY, envelope.MaxY) };
    }


    public double[] CentreOfEnveloppe(Geometry geom)
    {
        Envelope env = new Envelope();
        geom.GetEnvelope(env);
        return Centre(env);
    }

}
</code>







--
View this message in context: http://osgeo-org.1560.n6.nabble.com/gdal-dev-Ogr-Is-an-OGR-Geometry-InteriorPoint-method-tp4995621p4995701.html
Sent from the GDAL - Dev mailing list archive at Nabble.com.