Skip to content
Permalink
Browse files
PROJ6: more uses of reprojectionObj to improve performance
  • Loading branch information
rouault committed Oct 1, 2019
1 parent b73b82b commit 1ca010b949072d6de733d79fe3cc7e0e1bbdb28e
Show file tree
Hide file tree
Showing 12 changed files with 168 additions and 43 deletions.
@@ -201,7 +201,20 @@ int msComputeKernelDensityDataset(mapObj *map, imageObj *image, layerObj *kernel
values = (float*) msSmallCalloc(im_width * im_height, sizeof(float));
#ifdef USE_PROJ
if(layer->project)
msProjectShape(&layer->projection, &map->projection, &shape);
{
if( layer->reprojectorLayerToMap == NULL )
{
layer->reprojectorLayerToMap = msProjectCreateReprojector(
&layer->projection, &map->projection);
if( layer->reprojectorLayerToMap == NULL )
{
msFreeShape(&shape);
msLayerClose(layer);
return MS_FAILURE;
}
}
msProjectShapeEx(layer->reprojectorLayerToMap, &shape);
}
#endif

/* the weight for the sample is set to 1.0 by default. If the
@@ -251,7 +251,18 @@ int getNextShape(mapObj *map, layerObj *layer, double *values, int *nvalues, sty
if(status == MS_SUCCESS) {
#ifdef USE_PROJ
if(layer->project)
msProjectShape(&layer->projection, &map->projection, shape);
{
if( layer->reprojectorLayerToMap == NULL )
{
layer->reprojectorLayerToMap = msProjectCreateReprojector(
&layer->projection, &map->projection);
if( layer->reprojectorLayerToMap == NULL )
{
return MS_FAILURE;
}
}
msProjectShapeEx(layer->reprojectorLayerToMap, shape);
}
#endif

if(msBindLayerToShape(layer, shape, MS_DRAWMODE_FEATURES|MS_DRAWMODE_LABELS) != MS_SUCCESS)
@@ -1001,6 +1001,7 @@ int RebuildClusters(layerObj *layer, int isQuery)
#ifdef USE_CLUSTER_EXTERNAL
int layerIndex;
#endif
reprojectionObj* reprojector = NULL;

msClusterLayerInfo* layerinfo = layer->layerinfo;

@@ -1133,11 +1134,18 @@ int RebuildClusters(layerObj *layer, int isQuery)
if ((current = clusterInfoCreate(layerinfo)) == NULL)
return MS_FAILURE;

#if defined(USE_PROJ) && defined(USE_CLUSTER_EXTERNAL)
if(srcLayer->transform == MS_TRUE && srcLayer->project && layer->transform == MS_TRUE && layer->project &&msProjectionsDiffer(&(srcLayer->projection), &(layer->projection)))
{
reprojector = msProjectCreateReprojector(&srcLayer->projection, &layer->projection);
}
#endif

while((status = msLayerNextShape(srcLayer, &current->shape)) == MS_SUCCESS) {
#if defined(USE_PROJ) && defined(USE_CLUSTER_EXTERNAL)
/* transform the shape to the projection of this layer */
if(srcLayer->transform == MS_TRUE && srcLayer->project && layer->transform == MS_TRUE && layer->project &&msProjectionsDiffer(&(srcLayer->projection), &(layer->projection)))
msProjectShape(&srcLayer->projection, &layer->projection, &current->shape);
if( reprojector )
msProjectShapeEx(reprojector, &current->shape);
#endif
/* set up positions and variance */
current->avgx = current->x = current->shape.bounds.minx;
@@ -1173,6 +1181,7 @@ int RebuildClusters(layerObj *layer, int isQuery)
/* add this shape to the tree */
if (treeNodeAddShape(layerinfo, layerinfo->root, current, depth) != MS_SUCCESS) {
clusterInfoDestroyList(layerinfo, current);
msProjectDestroyReprojector(reprojector);
return MS_FAILURE;
}
}
@@ -1190,17 +1199,21 @@ int RebuildClusters(layerObj *layer, int isQuery)
/* if not found add this shape as a new cluster */
if (treeNodeAddShape(layerinfo, layerinfo->root, current, depth) != MS_SUCCESS) {
clusterInfoDestroyList(layerinfo, current);
msProjectDestroyReprojector(reprojector);
return MS_FAILURE;
}
}
}

if ((current = clusterInfoCreate(layerinfo)) == NULL) {
clusterInfoDestroyList(layerinfo, current);
msProjectDestroyReprojector(reprojector);
return MS_FAILURE;
}
}

msProjectDestroyReprojector(reprojector);

clusterInfoDestroyList(layerinfo, current);

if (layerinfo->algorithm == MSCLUSTER_ALGORITHM_FULL) {
@@ -1679,6 +1679,7 @@ int circleLayerDrawShape(mapObj *map, imageObj *image, layerObj *layer, shapeObj
/* TODO: need to handle circle annotation */
}

static
int pointLayerDrawShape(mapObj *map, imageObj *image, layerObj *layer, shapeObj *shape, int drawmode)
{
int l, c = shape->classindex, j, i, s;
@@ -1687,7 +1688,18 @@ int pointLayerDrawShape(mapObj *map, imageObj *image, layerObj *layer, shapeObj

#ifdef USE_PROJ
if (layer->project && layer->transform == MS_TRUE)
msProjectShape(&layer->projection, &map->projection, shape);
{
if( layer->reprojectorLayerToMap == NULL )
{
layer->reprojectorLayerToMap = msProjectCreateReprojector(
&layer->projection, &map->projection);
if( layer->reprojectorLayerToMap == NULL )
{
return MS_FAILURE;
}
}
msProjectShapeEx(layer->reprojectorLayerToMap, shape);
}
#endif

// Only take into account map rotation if the label and style angles are
@@ -2015,7 +2027,18 @@ int msDrawShape(mapObj *map, layerObj *layer, shapeObj *shape, imageObj *image,

#ifdef USE_PROJ
if (layer->project && layer->transform == MS_TRUE)
msProjectShape(&layer->projection, &map->projection, shape);
{
if( layer->reprojectorLayerToMap == NULL )
{
layer->reprojectorLayerToMap = msProjectCreateReprojector(
&layer->projection, &map->projection);
if( layer->reprojectorLayerToMap == NULL )
{
return MS_FAILURE;
}
}
msProjectShapeEx(layer->reprojectorLayerToMap, shape);
}
#endif

/* check if we'll need the unclipped shape */
@@ -3577,6 +3577,8 @@ int initLayer(layerObj *layer, mapObj *map)
}

layer->project = MS_TRUE;
layer->reprojectorLayerToMap = NULL;
layer->reprojectorMapToLayer = NULL;

initCluster(&layer->cluster);

@@ -3720,6 +3722,8 @@ int freeLayer(layerObj *layer)
msFree(layer->vtable);
msFree(layer->classgroup);

msProjectDestroyReprojector(layer->reprojectorLayerToMap);
msProjectDestroyReprojector(layer->reprojectorMapToLayer);
msFreeProjection(&(layer->projection));
msFreeExpression(&layer->_geomtransform);

@@ -702,8 +702,15 @@ graticuleIntersectionObj *msGraticuleLayerGetIntersectionPoints(mapObj *map,
msCopyShape(&shapegrid, &tmpshape);
/* status = msDrawShape(map, layer, &tmpshape, image, -1); */

if(layer->project)
msProjectShape(&layer->projection, &map->projection, &shapegrid);
if(layer->project) {
if( layer->reprojectorLayerToMap == NULL )
{
layer->reprojectorLayerToMap = msProjectCreateReprojector(
&layer->projection, &map->projection);
}
if( layer->reprojectorLayerToMap )
msProjectShapeEx(layer->reprojectorLayerToMap, &shapegrid);
}

msClipPolylineRect(&shapegrid, cliprect);

@@ -1070,7 +1077,13 @@ static int _AdjustLabelPosition( layerObj *pLayer, shapeObj *pShape, msGraticule
#ifdef USE_PROJ
if(pLayer->project)
{
msProjectShape( &pLayer->projection, &pLayer->map->projection, pShape );
if( pLayer->reprojectorLayerToMap == NULL )
{
pLayer->reprojectorLayerToMap = msProjectCreateReprojector(
&pLayer->projection, &pLayer->map->projection);
}
if( pLayer->reprojectorLayerToMap )
msProjectShapeEx(pLayer->reprojectorLayerToMap, pShape );

/* Poor man detection of reprojection failure */
if( msProjIsGeographicCRS(&(pLayer->projection)) !=
@@ -1191,7 +1204,13 @@ static int _AdjustLabelPosition( layerObj *pLayer, shapeObj *pShape, msGraticule
}
}

msProjectShape( &pLayer->map->projection, &pLayer->projection, pShape );
if( pLayer->reprojectorMapToLayer == NULL )
{
pLayer->reprojectorMapToLayer = msProjectCreateReprojector(
&pLayer->map->projection, &pLayer->projection);
}
if( pLayer->reprojectorMapToLayer )
msProjectShapeEx(pLayer->reprojectorMapToLayer, pShape );
}
#endif

@@ -1589,6 +1589,7 @@ int LayerDefaultGetShapeCount(layerObj *layer, rectObj rect, projectionObj *rect
shapeObj shape, searchshape;
int nShapeCount = 0;
rectObj searchrect = rect;
reprojectionObj* reprojector = NULL;

msInitShape(&searchshape);
msRectToPolygon(searchrect, &searchshape);
@@ -1622,7 +1623,12 @@ int LayerDefaultGetShapeCount(layerObj *layer, rectObj rect, projectionObj *rect
{
#ifdef USE_PROJ
if(layer->project && msProjectionsDiffer(&(layer->projection), rectProjection))
msProjectShape(&(layer->projection), rectProjection, &shape);
{
if( reprojector == NULL )
reprojector = msProjectCreateReprojector(&(layer->projection), rectProjection);
if( reprojector )
msProjectShapeEx(reprojector, &shape);
}
else
layer->project = MS_FALSE;
#endif
@@ -1656,6 +1662,7 @@ int LayerDefaultGetShapeCount(layerObj *layer, rectObj rect, projectionObj *rect
}

msFreeShape(&searchshape);
msProjectDestroyReprojector(reprojector);

return nShapeCount;
}
@@ -510,8 +510,17 @@ int msMVTWriteTile( mapObj *map, int sendheaders ) {
}

if( layer->project ) {
status = msProjectShape(&layer->projection, &layer->map->projection, &shape);
} if( status == MS_SUCCESS ) {
if( layer->reprojectorLayerToMap == NULL )
{
layer->reprojectorLayerToMap = msProjectCreateReprojector(
&layer->projection, &map->projection);
}
if( layer->reprojectorLayerToMap )
status = msProjectShapeEx(layer->reprojectorLayerToMap, &shape);
else
status = MS_FAILURE;
}
if( status == MS_SUCCESS ) {
status = mvtWriteShape( layer, &shape, mvt_layer, item_list, &value_lookup_cache, &map->extent, buffer );
}

@@ -1179,9 +1179,15 @@ int msOGRWriteFromQuery( mapObj *map, outputFormatObj *format, int sendheaders )
}

if( layer->project ) {
status =
msProjectShape(&layer->projection, &layer->map->projection,
&resultshape);
if( layer->reprojectorLayerToMap == NULL )
{
layer->reprojectorLayerToMap = msProjectCreateReprojector(
&layer->projection, &layer->map->projection);
}
if( layer->reprojectorLayerToMap )
status = msProjectShapeEx(layer->reprojectorLayerToMap, &resultshape);
else
status = MS_FAILURE;
}

/*
@@ -40,8 +40,7 @@
static char *ms_proj_lib = NULL;

static int msTestNeedWrap( pointObj pt1, pointObj pt2, pointObj pt2_geo,
projectionObj *src_proj,
projectionObj *dst_proj );
reprojectionObj* reprojector );

#if defined(USE_PROJ) && PROJ_VERSION_MAJOR >= 6

@@ -964,7 +963,7 @@ static void msProjectGrowRect(reprojectionObj* reprojector,
/* reprojects and the other end does not. Finds the horizon. */
/************************************************************************/
#ifdef USE_PROJ
static int msProjectSegment( projectionObj *in, projectionObj *out,
static int msProjectSegment( reprojectionObj* reprojector,
pointObj *start, pointObj *end )

{
@@ -976,12 +975,12 @@ static int msProjectSegment( projectionObj *in, projectionObj *out,
/* then re-call with the points reversed. */
/* -------------------------------------------------------------------- */
testPoint = *start;
if( msProjectPoint( in, out, &testPoint ) == MS_FAILURE ) {
if( msProjectPointEx( reprojector, &testPoint ) == MS_FAILURE ) {
testPoint = *end;
if( msProjectPoint( in, out, &testPoint ) == MS_FAILURE )
if( msProjectPointEx( reprojector, &testPoint ) == MS_FAILURE )
return MS_FAILURE;
else
return msProjectSegment( in, out, end, start );
return msProjectSegment( reprojector, end, start );
}

/* -------------------------------------------------------------------- */
@@ -1002,7 +1001,7 @@ static int msProjectSegment( projectionObj *in, projectionObj *out,

testPoint = midPoint;

if( msProjectPoint( in, out, &testPoint ) == MS_FAILURE ) {
if( msProjectPointEx( reprojector, &testPoint ) == MS_FAILURE ) {
if (midPoint.x == subEnd.x && midPoint.y == subEnd.y)
return MS_FAILURE; /* break infinite loop */

@@ -1021,8 +1020,8 @@ static int msProjectSegment( projectionObj *in, projectionObj *out,
/* -------------------------------------------------------------------- */
*end = subStart;

if( msProjectPoint( in, out, end ) == MS_FAILURE
|| msProjectPoint( in, out, start ) == MS_FAILURE )
if( msProjectPointEx( reprojector, end ) == MS_FAILURE
|| msProjectPointEx( reprojector, start ) == MS_FAILURE )
return MS_FAILURE;
else
return MS_SUCCESS;
@@ -1126,7 +1125,7 @@ msProjectShapeLine(reprojectionObj* reprojector,
dist = wrkPoint.x - pt1Geo.x;
if( fabs(dist) > 180.0
&& msTestNeedWrap( thisPoint, lastPoint,
pt1Geo, in, out ) ) {
pt1Geo, reprojector ) ) {
if( dist > 0.0 )
wrkPoint.x -= 360.0;
else if( dist < 0.0 )
@@ -1150,7 +1149,7 @@ msProjectShapeLine(reprojectionObj* reprojector,

startPoint = lastPoint;
endPoint = thisPoint;
if( msProjectSegment( in, out, &startPoint, &endPoint )
if( msProjectSegment( reprojector, &startPoint, &endPoint )
== MS_SUCCESS ) {
line_out->point[line_out->numpoints++] = endPoint;
}
@@ -1181,7 +1180,7 @@ msProjectShapeLine(reprojectionObj* reprojector,

startPoint = lastPoint;
endPoint = thisPoint;
if( msProjectSegment( in, out, &endPoint, &startPoint )
if( msProjectSegment( reprojector, &endPoint, &startPoint )
== MS_SUCCESS ) {
lineObj newLine;

@@ -1382,7 +1381,7 @@ int msProjectLineEx(reprojectionObj* reprojector, lineObj *line)
dist = line->point[i].x - line->point[0].x;
if( fabs(dist) > 180.0 ) {
if( msTestNeedWrap( thisPoint, startPoint,
line->point[0], reprojector->in, reprojector->out ) ) {
line->point[0], reprojector ) ) {
if( dist > 0.0 ) {
line->point[i].x -= 360.0;
} else if( dist < 0.0 ) {
@@ -2125,18 +2124,17 @@ reprojected to geographic space. However, it does not:

#ifdef USE_PROJ
static int msTestNeedWrap( pointObj pt1, pointObj pt2, pointObj pt2_geo,
projectionObj *in,
projectionObj *out )
reprojectionObj* reprojector )

{
pointObj middle;

middle.x = (pt1.x + pt2.x) * 0.5;
middle.y = (pt1.y + pt2.y) * 0.5;

if( msProjectPoint( in, out, &pt1 ) == MS_FAILURE
|| msProjectPoint( in, out, &pt2 ) == MS_FAILURE
|| msProjectPoint( in, out, &middle ) == MS_FAILURE )
if( msProjectPointEx( reprojector, &pt1 ) == MS_FAILURE
|| msProjectPointEx( reprojector, &pt2 ) == MS_FAILURE
|| msProjectPointEx( reprojector, &middle ) == MS_FAILURE )
return 0;

/*

0 comments on commit 1ca010b

Please sign in to comment.