www.pudn.com > mitab.rar > mitab_geometry.cpp
/**********************************************************************
* $Id: mitab_geometry.cpp,v 1.5 2004/06/30 20:29:04 dmorissette Exp $
*
* Name: mitab_geometry.cpp
* Project: MapInfo TAB Read/Write library
* Language: C++
* Purpose: Geometry manipulation functions.
* Author: Daniel Morissette, dmorissette@dmsolutions.ca
* Based on functions from mapprimitive.c/mapsearch.c in the source
* of UMN MapServer by Stephen Lime (http://mapserver.gis.umn.edu/)
**********************************************************************
* Copyright (c) 1999-2001, Daniel Morissette
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
**********************************************************************
*
* $Log: mitab_geometry.cpp,v $
* Revision 1.5 2004/06/30 20:29:04 dmorissette
* Fixed refs to old address danmo@videotron.ca
*
* Revision 1.4 2001/12/18 23:42:28 daniel
* Added a test in OGRPolygonLabelPoint() to prevent returning a point
* outside of the polygon MBR (bug 673).
*
* Revision 1.3 2001/01/22 16:03:58 warmerda
* expanded tabs
*
* Revision 1.2 2000/09/28 16:39:44 warmerda
* avoid warnings for unused, and unitialized variables
*
* Revision 1.1 2000/09/19 17:19:40 daniel
* Initial Revision
*
**********************************************************************/
#include "ogr_geometry.h"
#define OGR_NUM_RINGS(poly) (poly->getNumInteriorRings()+1)
#define OGR_GET_RING(poly, i) (i==0?poly->getExteriorRing():poly->getInteriorRing(i-1))
/**********************************************************************
* OGRPointInRing()
*
* Returns TRUE is point is inside ring, FALSE otherwise
*
* Adapted version of msPointInPolygon() from MapServer's mapsearch.c
**********************************************************************/
GBool OGRPointInRing(OGRPoint *poPoint, OGRLineString *poRing)
{
int i, j, numpoints;
GBool status = FALSE;
double x, y;
numpoints = poRing->getNumPoints();
x = poPoint->getX();
y = poPoint->getY();
for (i = 0, j = numpoints-1; i < numpoints; j = i++)
{
if ((((poRing->getY(i)<=y) && (ygetY(j))) ||
((poRing->getY(j)<=y) && (ygetY(i)))) &&
(x < (poRing->getX(j) - poRing->getX(i)) * (y - poRing->getY(i)) /
(poRing->getY(j) - poRing->getY(i)) + poRing->getX(i)))
status = !status;
}
return status;
}
/**********************************************************************
* OGRIntersectPointPolygon()
*
* Instead of using ring orientation we count the number of parts the
* point falls in. If odd the point is in the polygon, if 0 or even
* then the point is in a hole or completely outside.
*
* Returns TRUE is point is inside polygon, FALSE otherwise
*
* Adapted version of msIntersectPointPolygon() from MapServer's mapsearch.c
**********************************************************************/
GBool OGRIntersectPointPolygon(OGRPoint *poPoint, OGRPolygon *poPoly)
{
int i;
GBool status = FALSE;
for(i=0; i (max) ? CLIP_RIGHT : CLIP_MIDDLE));
#define ROUND(a) ( (a) + 0.5 )
#define SWAP( a, b, t) ( (t) = (a), (a) = (b), (b) = (t) )
#define EDGE_CHECK( x0, x, x1) ((x) < MIN( (x0), (x1)) ? CLIP_LEFT : ((x) > MAX( (x0), (x1)) ? CLIP_RIGHT : CLIP_MIDDLE ))
#define NUM_SCANLINES 5
int OGRPolygonLabelPoint(OGRPolygon *poPoly, OGRPoint *poLabelPoint)
{
double slope;
OGRRawPoint point1, point2;
int i, j, k, nfound;
double x, y, *xintersect, temp;
double hi_y, lo_y;
int wrong_order, n;
double len, max_len=0;
double skip;
OGREnvelope oEnv;
if (poPoly == NULL)
return OGRERR_FAILURE;
poPoly->getEnvelope(&oEnv);
poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);
poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);
//if(get_centroid(p, lp, &miny, &maxy) == -1) return(-1);
if(OGRIntersectPointPolygon(poLabelPoint, poPoly) == TRUE) /* cool, done */
return OGRERR_NONE;
/* do it the hard way - scanline */
skip = (oEnv.MaxY - oEnv.MinY)/NUM_SCANLINES;
n=0;
for(j=0; jgetNumPoints();
}
xintersect = (double *)calloc(n, sizeof(double));
if (xintersect == NULL)
return OGRERR_FAILURE;
for(k=1; k<=NUM_SCANLINES; k++)
{
/* sample the shape in the y direction */
y = oEnv.MaxY - k*skip;
/* need to find a y that won't intersect any vertices exactly */
hi_y = y - 1; /* first initializing lo_y, hi_y to be any 2 pnts on either side of lp->y */
lo_y = y + 1;
for(j=0; j= y))
break; /* already initialized */
for(i=0; i < poRing->getNumPoints(); i++)
{
if((lo_y < y) && (hi_y >= y))
break; /* already initialized */
if(poRing->getY(i) < y)
lo_y = poRing->getY(i);
if(poRing->getY(i) >= y)
hi_y = poRing->getY(i);
}
}
n=0;
for(j=0; jgetNumPoints(); i++)
{
if((poRing->getY(i) < y) &&
((y - poRing->getY(i)) < (y - lo_y)))
lo_y = poRing->getY(i);
if((poRing->getY(i) >= y) &&
((poRing->getY(i) - y) < (hi_y - y)))
hi_y = poRing->getY(i);
}
}
if(lo_y == hi_y)
return OGRERR_FAILURE;
else
y = (hi_y + lo_y)/2.0;
nfound = 0;
for(j=0; jgetX(poRing->getNumPoints()-1);
point1.y = poRing->getY(poRing->getNumPoints()-1);
for(i=0; i < poRing->getNumPoints(); i++)
{
point2.x = poRing->getX(i);
point2.y = poRing->getY(i);
if(EDGE_CHECK(point1.y, y, point2.y) == CLIP_MIDDLE)
{
if(point1.y == point2.y)
continue; /* ignore horizontal edges */
else
slope = (point2.x - point1.x) / (point2.y - point1.y);
x = point1.x + (y - point1.y)*slope;
xintersect[nfound++] = x;
} /* End of checking this edge */
point1 = point2; /* Go on to next edge */
}
} /* Finished the scanline */
/* First, sort the intersections */
do
{
wrong_order = 0;
for(i=0; i < nfound-1; i++)
{
if(xintersect[i] > xintersect[i+1])
{
wrong_order = 1;
SWAP(xintersect[i], xintersect[i+1], temp);
}
}
} while(wrong_order);
/* Great, now find longest span */
point1.y = point2.y = y;
for(i=0; i < nfound; i += 2)
{
point1.x = xintersect[i];
point2.x = xintersect[i+1];
/* len = length(point1, point2); */
len = ABS((point2.x - point1.x));
if(len > max_len)
{
max_len = len;
poLabelPoint->setX( (point1.x + point2.x)/2 );
poLabelPoint->setY( y );
}
}
}
free(xintersect);
/* __TODO__ Bug 673
* There seem to be some polygons for which the label is returned
* completely outside of the polygon's MBR and this messes the
* file bounds, etc.
* Until we find the source of the problem, we'll at least validate
* the label point to make sure that it overlaps the polygon MBR.
*/
if( poLabelPoint->getX() < oEnv.MinX
|| poLabelPoint->getY() < oEnv.MinY
|| poLabelPoint->getX() > oEnv.MaxX
|| poLabelPoint->getY() > oEnv.MaxY )
{
// Reset label coordinates to center of MBR, just in case
poLabelPoint->setX((oEnv.MaxX + oEnv.MinX)/2.0);
poLabelPoint->setY((oEnv.MaxY + oEnv.MinY)/2.0);
// And return an error
return OGRERR_FAILURE;
}
if(max_len > 0)
return OGRERR_NONE;
else
return OGRERR_FAILURE;
}
/**********************************************************************
* OGRGetCentroid()
*
* Calculate polygon gravity center.
*
* Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise.
*
* Adapted version of get_centroid() from MapServer's mapprimitive.c
**********************************************************************/
int OGRGetCentroid(OGRPolygon *poPoly, OGRPoint *poCentroid)
{
int i,j;
double cent_weight_x=0.0, cent_weight_y=0.0;
double len, total_len=0;
for(i=0; igetX(0);
y2 = poRing->getY(0);
for(j=1; jgetNumPoints(); j++)
{
x1 = x2;
y1 = y2;
x2 = poRing->getX(j);
y2 = poRing->getY(j);
len = sqrt( pow((x2-x1),2) + pow((y2-y1),2) );
cent_weight_x += len * ((x1 + x2)/2.0);
cent_weight_y += len * ((y1 + y2)/2.0);
total_len += len;
}
}
if(total_len == 0)
return(OGRERR_FAILURE);
poCentroid->setX( cent_weight_x / total_len );
poCentroid->setY( cent_weight_y / total_len );
return OGRERR_NONE;
}
/**********************************************************************
* OGRPolylineCenterPoint()
*
* Return the center point of a polyline.
*
* In MapInfo, for a simple or multiple polyline (pline), the center point
* in the object definition is supposed to be either the center point of
* the pline or the first section of a multiple pline (if an odd number of
* points in the pline or first section), or the midway point between the
* two central points (if an even number of points involved).
*
* Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise.
**********************************************************************/
int OGRPolylineCenterPoint(OGRLineString *poLine, OGRPoint *poLabelPoint)
{
if (poLine == NULL || poLine->getNumPoints() < 2)
return OGRERR_FAILURE;
if (poLine->getNumPoints() % 2 == 0)
{
// Return the midway between the 2 center points
int i = poLine->getNumPoints()/2;
poLabelPoint->setX( (poLine->getX(i-1) + poLine->getX(i))/2.0 );
poLabelPoint->setY( (poLine->getY(i-1) + poLine->getY(i))/2.0 );
}
else
{
// Return the center point
poLine->getPoint(poLine->getNumPoints()/2, poLabelPoint);
}
return OGRERR_NONE;
}
/**********************************************************************
* OGRPolylineLabelPoint()
*
* Generate a label point on a polyline: The center of the longest segment.
*
* Returns OGRERR_NONE if it succeeds or OGRERR_FAILURE otherwise.
**********************************************************************/
int OGRPolylineLabelPoint(OGRLineString *poLine, OGRPoint *poLabelPoint)
{
double segment_length, max_segment_length = 0.0;
double x1, y1, x2, y2;
if (poLine == NULL || poLine->getNumPoints() < 2)
return OGRERR_FAILURE;
max_segment_length = -1.0;
x2 = poLine->getX(0);
y2 = poLine->getY(0);
for(int i=1; igetNumPoints(); i++)
{
x1 = x2;
y1 = y2;
x2 = poLine->getX(i);
y2 = poLine->getY(i);
segment_length = pow((x2-x1),2) + pow((y2-y1),2);
if (segment_length > max_segment_length)
{
max_segment_length = segment_length;
poLabelPoint->setX( (x1 + x2)/2.0 );
poLabelPoint->setY( (y1 + y2)/2.0 );
}
}
return OGRERR_NONE;
}