Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  High Performance Image Reprojector
       4             :  * Purpose:  Implement cutline/blend mask generator.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2008, Frank Warmerdam <warmerdam@pobox.com>
       9             :  * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * Permission is hereby granted, free of charge, to any person obtaining a
      12             :  * copy of this software and associated documentation files (the "Software"),
      13             :  * to deal in the Software without restriction, including without limitation
      14             :  * the rights to use, copy, modify, merge, publish, distribute, sublicense,
      15             :  * and/or sell copies of the Software, and to permit persons to whom the
      16             :  * Software is furnished to do so, subject to the following conditions:
      17             :  *
      18             :  * The above copyright notice and this permission notice shall be included
      19             :  * in all copies or substantial portions of the Software.
      20             :  *
      21             :  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
      22             :  * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      23             :  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
      24             :  * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      25             :  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      26             :  * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      27             :  * DEALINGS IN THE SOFTWARE.
      28             :  ****************************************************************************/
      29             : 
      30             : #include "cpl_port.h"
      31             : #include "gdalwarper.h"
      32             : 
      33             : #include <cmath>
      34             : #include <cstdio>
      35             : #include <cstring>
      36             : #include <algorithm>
      37             : 
      38             : #include "cpl_conv.h"
      39             : #include "cpl_error.h"
      40             : #include "cpl_string.h"
      41             : #include "gdal.h"
      42             : #include "gdal_alg.h"
      43             : #include "gdal_priv.h"
      44             : #include "memdataset.h"
      45             : #include "ogr_api.h"
      46             : #include "ogr_core.h"
      47             : #include "ogr_geometry.h"
      48             : #include "ogr_geos.h"
      49             : 
      50             : /************************************************************************/
      51             : /*                         BlendMaskGenerator()                         */
      52             : /************************************************************************/
      53             : 
      54             : #ifndef HAVE_GEOS
      55             : 
      56             : static CPLErr BlendMaskGenerator(int /* nXOff */, int /* nYOff */,
      57             :                                  int /* nXSize */, int /* nYSize */,
      58             :                                  GByte * /* pabyPolyMask */,
      59             :                                  float * /* pafValidityMask */,
      60             :                                  OGRGeometryH /* hPolygon */,
      61             :                                  double /* dfBlendDist */)
      62             : {
      63             :     CPLError(CE_Failure, CPLE_AppDefined,
      64             :              "Blend distance support not available without the GEOS library.");
      65             :     return CE_Failure;
      66             : }
      67             : #else
      68           1 : static CPLErr BlendMaskGenerator(int nXOff, int nYOff, int nXSize, int nYSize,
      69             :                                  GByte *pabyPolyMask, float *pafValidityMask,
      70             :                                  OGRGeometryH hPolygon, double dfBlendDist)
      71             : {
      72             : 
      73             :     /* -------------------------------------------------------------------- */
      74             :     /*      Convert the polygon into a collection of lines so that we       */
      75             :     /*      measure distance from the edge even on the inside.              */
      76             :     /* -------------------------------------------------------------------- */
      77           1 :     const auto poPolygon = OGRGeometry::FromHandle(hPolygon);
      78             :     auto poLines = std::unique_ptr<OGRGeometry>(
      79           2 :         OGRGeometryFactory::forceToMultiLineString(poPolygon->clone()));
      80             : 
      81             :     /* -------------------------------------------------------------------- */
      82             :     /*      Prepare a clipping polygon a bit bigger than the area of        */
      83             :     /*      interest in the hopes of simplifying the cutline down to        */
      84             :     /*      stuff that will be relevant for this area of interest.          */
      85             :     /* -------------------------------------------------------------------- */
      86           2 :     CPLString osClipRectWKT;
      87             : 
      88             :     osClipRectWKT.Printf(
      89           1 :         "POLYGON((%g %g,%g %g,%g %g,%g %g,%g %g))", nXOff - (dfBlendDist + 1),
      90           1 :         nYOff - (dfBlendDist + 1), nXOff + nXSize + (dfBlendDist + 1),
      91           1 :         nYOff - (dfBlendDist + 1), nXOff + nXSize + (dfBlendDist + 1),
      92           1 :         nYOff + nYSize + (dfBlendDist + 1), nXOff - (dfBlendDist + 1),
      93           1 :         nYOff + nYSize + (dfBlendDist + 1), nXOff - (dfBlendDist + 1),
      94           1 :         nYOff - (dfBlendDist + 1));
      95             : 
      96           2 :     OGRPolygon oClipRect;
      97             :     {
      98           1 :         const char *pszClipRectWKT = osClipRectWKT.c_str();
      99           1 :         oClipRect.importFromWkt(&pszClipRectWKT);
     100             :     }
     101             : 
     102             :     // If it does not intersect the polym, zero the mask and return.
     103           1 :     if (!poPolygon->Intersects(&oClipRect))
     104             :     {
     105           0 :         memset(pafValidityMask, 0, sizeof(float) * nXSize * nYSize);
     106             : 
     107           0 :         return CE_None;
     108             :     }
     109             : 
     110             :     // If it does not intersect the line at all, just return.
     111           1 :     else if (!poLines->Intersects(&oClipRect))
     112             :     {
     113             : 
     114           0 :         return CE_None;
     115             :     }
     116             : 
     117           1 :     poLines.reset(poLines->Intersection(&oClipRect));
     118             : 
     119             :     /* -------------------------------------------------------------------- */
     120             :     /*      Convert our polygon into GEOS format, and compute an            */
     121             :     /*      envelope to accelerate later distance operations.               */
     122             :     /* -------------------------------------------------------------------- */
     123           1 :     OGREnvelope sEnvelope;
     124           1 :     GEOSContextHandle_t hGEOSCtxt = OGRGeometry::createGEOSContext();
     125           1 :     GEOSGeom poGEOSPoly = poLines->exportToGEOS(hGEOSCtxt);
     126           1 :     OGR_G_GetEnvelope(hPolygon, &sEnvelope);
     127             : 
     128             :     // This check was already done in the calling
     129             :     // function and should never be true.
     130             : 
     131             :     // if( sEnvelope.MinY - dfBlendDist > nYOff+nYSize
     132             :     //     || sEnvelope.MaxY + dfBlendDist < nYOff
     133             :     //     || sEnvelope.MinX - dfBlendDist > nXOff+nXSize
     134             :     //     || sEnvelope.MaxX + dfBlendDist < nXOff )
     135             :     //     return CE_None;
     136             : 
     137             :     const int iXMin = std::max(
     138           1 :         0, static_cast<int>(floor(sEnvelope.MinX - dfBlendDist - nXOff)));
     139             :     const int iXMax = std::min(
     140           1 :         nXSize, static_cast<int>(ceil(sEnvelope.MaxX + dfBlendDist - nXOff)));
     141             :     const int iYMin = std::max(
     142           1 :         0, static_cast<int>(floor(sEnvelope.MinY - dfBlendDist - nYOff)));
     143             :     const int iYMax = std::min(
     144           1 :         nYSize, static_cast<int>(ceil(sEnvelope.MaxY + dfBlendDist - nYOff)));
     145             : 
     146             :     /* -------------------------------------------------------------------- */
     147             :     /*      Loop over potential area within blend line distance,            */
     148             :     /*      processing each pixel.                                          */
     149             :     /* -------------------------------------------------------------------- */
     150         101 :     for (int iY = 0; iY < nYSize; iY++)
     151             :     {
     152         100 :         double dfLastDist = 0.0;
     153             : 
     154       10100 :         for (int iX = 0; iX < nXSize; iX++)
     155             :         {
     156       10000 :             if (iX < iXMin || iX >= iXMax || iY < iYMin || iY > iYMax ||
     157        3060 :                 dfLastDist > dfBlendDist + 1.5)
     158             :             {
     159        7944 :                 if (pabyPolyMask[iX + iY * nXSize] == 0)
     160        7784 :                     pafValidityMask[iX + iY * nXSize] = 0.0;
     161             : 
     162        7944 :                 dfLastDist -= 1.0;
     163        8528 :                 continue;
     164             :             }
     165             : 
     166        2056 :             CPLString osPointWKT;
     167        2056 :             osPointWKT.Printf("POINT(%d.5 %d.5)", iX + nXOff, iY + nYOff);
     168             : 
     169        2056 :             GEOSGeom poGEOSPoint = GEOSGeomFromWKT_r(hGEOSCtxt, osPointWKT);
     170             : 
     171        2056 :             double dfDist = 0.0;
     172        2056 :             GEOSDistance_r(hGEOSCtxt, poGEOSPoly, poGEOSPoint, &dfDist);
     173        2056 :             GEOSGeom_destroy_r(hGEOSCtxt, poGEOSPoint);
     174             : 
     175        2056 :             dfLastDist = dfDist;
     176             : 
     177        2056 :             if (dfDist > dfBlendDist)
     178             :             {
     179         584 :                 if (pabyPolyMask[iX + iY * nXSize] == 0)
     180         366 :                     pafValidityMask[iX + iY * nXSize] = 0.0;
     181             : 
     182         584 :                 continue;
     183             :             }
     184             : 
     185        1472 :             const double dfRatio =
     186        1472 :                 pabyPolyMask[iX + iY * nXSize] == 0
     187        1472 :                     ? 0.5 - (dfDist / dfBlendDist) * 0.5   // Outside.
     188         622 :                     : 0.5 + (dfDist / dfBlendDist) * 0.5;  // Inside.
     189             : 
     190        1472 :             pafValidityMask[iX + iY * nXSize] *= static_cast<float>(dfRatio);
     191             :         }
     192             :     }
     193             : 
     194             :     /* -------------------------------------------------------------------- */
     195             :     /*      Cleanup                                                         */
     196             :     /* -------------------------------------------------------------------- */
     197           1 :     GEOSGeom_destroy_r(hGEOSCtxt, poGEOSPoly);
     198           1 :     OGRGeometry::freeGEOSContext(hGEOSCtxt);
     199             : 
     200           1 :     return CE_None;
     201             : }
     202             : #endif  // HAVE_GEOS
     203             : 
     204             : /************************************************************************/
     205             : /*                         CutlineTransformer()                         */
     206             : /*                                                                      */
     207             : /*      A simple transformer for the cutline that just offsets          */
     208             : /*      relative to the current chunk.                                  */
     209             : /************************************************************************/
     210             : 
     211          29 : static int CutlineTransformer(void *pTransformArg, int bDstToSrc,
     212             :                               int nPointCount, double *x, double *y,
     213             :                               double * /* z */, int * /* panSuccess */)
     214             : {
     215          29 :     int nXOff = static_cast<int *>(pTransformArg)[0];
     216          29 :     int nYOff = static_cast<int *>(pTransformArg)[1];
     217             : 
     218          29 :     if (bDstToSrc)
     219             :     {
     220           0 :         nXOff *= -1;
     221           0 :         nYOff *= -1;
     222             :     }
     223             : 
     224        5005 :     for (int i = 0; i < nPointCount; i++)
     225             :     {
     226        4976 :         x[i] -= nXOff;
     227        4976 :         y[i] -= nYOff;
     228             :     }
     229             : 
     230          29 :     return TRUE;
     231             : }
     232             : 
     233             : /************************************************************************/
     234             : /*                       GDALWarpCutlineMasker()                        */
     235             : /*                                                                      */
     236             : /*      This function will generate a source mask based on a            */
     237             : /*      provided cutline, and optional blend distance.                  */
     238             : /************************************************************************/
     239             : 
     240           0 : CPLErr GDALWarpCutlineMasker(void *pMaskFuncArg, int nBandCount,
     241             :                              GDALDataType eType, int nXOff, int nYOff,
     242             :                              int nXSize, int nYSize, GByte **ppImageData,
     243             :                              int bMaskIsFloat, void *pValidityMask)
     244             : 
     245             : {
     246           0 :     return GDALWarpCutlineMaskerEx(pMaskFuncArg, nBandCount, eType, nXOff,
     247             :                                    nYOff, nXSize, nYSize, ppImageData,
     248           0 :                                    bMaskIsFloat, pValidityMask, nullptr);
     249             : }
     250             : 
     251          33 : CPLErr GDALWarpCutlineMaskerEx(void *pMaskFuncArg, int /* nBandCount */,
     252             :                                GDALDataType /* eType */, int nXOff, int nYOff,
     253             :                                int nXSize, int nYSize,
     254             :                                GByte ** /*ppImageData */, int bMaskIsFloat,
     255             :                                void *pValidityMask, int *pnValidityFlag)
     256             : 
     257             : {
     258          33 :     if (pnValidityFlag)
     259          33 :         *pnValidityFlag = GCMVF_PARTIAL_INTERSECTION;
     260             : 
     261          33 :     if (nXSize < 1 || nYSize < 1)
     262           0 :         return CE_None;
     263             : 
     264             :     /* -------------------------------------------------------------------- */
     265             :     /*      Do some minimal checking.                                       */
     266             :     /* -------------------------------------------------------------------- */
     267          33 :     if (!bMaskIsFloat)
     268             :     {
     269           0 :         CPLAssert(false);
     270             :         return CE_Failure;
     271             :     }
     272             : 
     273          33 :     GDALWarpOptions *psWO = static_cast<GDALWarpOptions *>(pMaskFuncArg);
     274             : 
     275          33 :     if (psWO == nullptr || psWO->hCutline == nullptr)
     276             :     {
     277           0 :         CPLAssert(false);
     278             :         return CE_Failure;
     279             :     }
     280             : 
     281          33 :     GDALDriverH hMemDriver = GDALGetDriverByName("MEM");
     282          33 :     if (hMemDriver == nullptr)
     283             :     {
     284           0 :         CPLError(CE_Failure, CPLE_AppDefined,
     285             :                  "GDALWarpCutlineMasker needs MEM driver");
     286           0 :         return CE_Failure;
     287             :     }
     288             : 
     289             :     /* -------------------------------------------------------------------- */
     290             :     /*      Check the polygon.                                              */
     291             :     /* -------------------------------------------------------------------- */
     292          33 :     OGRGeometryH hPolygon = static_cast<OGRGeometryH>(psWO->hCutline);
     293             : 
     294          63 :     if (wkbFlatten(OGR_G_GetGeometryType(hPolygon)) != wkbPolygon &&
     295          30 :         wkbFlatten(OGR_G_GetGeometryType(hPolygon)) != wkbMultiPolygon)
     296             :     {
     297           0 :         CPLError(CE_Failure, CPLE_NotSupported,
     298             :                  "Cutline should be a polygon or a multipolygon");
     299           0 :         return CE_Failure;
     300             :     }
     301             : 
     302          33 :     OGREnvelope sEnvelope;
     303          33 :     OGR_G_GetEnvelope(hPolygon, &sEnvelope);
     304             : 
     305          33 :     float *pafMask = static_cast<float *>(pValidityMask);
     306             : 
     307          33 :     if (sEnvelope.MaxX + psWO->dfCutlineBlendDist < nXOff ||
     308          33 :         sEnvelope.MinX - psWO->dfCutlineBlendDist > nXOff + nXSize ||
     309          33 :         sEnvelope.MaxY + psWO->dfCutlineBlendDist < nYOff ||
     310          33 :         sEnvelope.MinY - psWO->dfCutlineBlendDist > nYOff + nYSize)
     311             :     {
     312           0 :         if (pnValidityFlag)
     313           0 :             *pnValidityFlag = GCMVF_NO_INTERSECTION;
     314             : 
     315             :         // We are far from the blend line - everything is masked to zero.
     316             :         // It would be nice to realize no work is required for this whole
     317             :         // chunk!
     318           0 :         memset(pafMask, 0, sizeof(float) * nXSize * nYSize);
     319           0 :         return CE_None;
     320             :     }
     321             : 
     322             :     // And now check if the chunk to warp is fully contained within the cutline
     323             :     // to save rasterization.
     324          33 :     if (OGRGeometryFactory::haveGEOS()
     325             : #ifdef DEBUG
     326             :         // Env var just for debugging purposes
     327          33 :         && !CPLTestBool(
     328             :                CPLGetConfigOption("GDALCUTLINE_SKIP_CONTAINMENT_TEST", "NO"))
     329             : #endif
     330             :     )
     331             :     {
     332          33 :         OGRLinearRing *poRing = new OGRLinearRing();
     333          33 :         poRing->addPoint(-psWO->dfCutlineBlendDist + nXOff,
     334          33 :                          -psWO->dfCutlineBlendDist + nYOff);
     335          33 :         poRing->addPoint(-psWO->dfCutlineBlendDist + nXOff,
     336          33 :                          psWO->dfCutlineBlendDist + nYOff + nYSize);
     337          33 :         poRing->addPoint(psWO->dfCutlineBlendDist + nXOff + nXSize,
     338          33 :                          psWO->dfCutlineBlendDist + nYOff + nYSize);
     339          33 :         poRing->addPoint(psWO->dfCutlineBlendDist + nXOff + nXSize,
     340          33 :                          -psWO->dfCutlineBlendDist + nYOff);
     341          33 :         poRing->addPoint(-psWO->dfCutlineBlendDist + nXOff,
     342          33 :                          -psWO->dfCutlineBlendDist + nYOff);
     343          33 :         OGRPolygon oChunkFootprint;
     344          33 :         oChunkFootprint.addRingDirectly(poRing);
     345          33 :         OGREnvelope sChunkEnvelope;
     346          33 :         oChunkFootprint.getEnvelope(&sChunkEnvelope);
     347          42 :         if (sEnvelope.Contains(sChunkEnvelope) &&
     348           9 :             OGRGeometry::FromHandle(hPolygon)->Contains(&oChunkFootprint))
     349             :         {
     350           8 :             if (pnValidityFlag)
     351           8 :                 *pnValidityFlag = GCMVF_CHUNK_FULLY_WITHIN_CUTLINE;
     352             : 
     353           8 :             CPLDebug("WARP", "Source chunk fully contained within cutline.");
     354           8 :             return CE_None;
     355             :         }
     356             :     }
     357             : 
     358             :     /* -------------------------------------------------------------------- */
     359             :     /*      Create a byte buffer into which we can burn the                 */
     360             :     /*      mask polygon and wrap it up as a memory dataset.                */
     361             :     /* -------------------------------------------------------------------- */
     362          25 :     GByte *pabyPolyMask = static_cast<GByte *>(CPLCalloc(nXSize, nYSize));
     363             : 
     364             :     auto poMEMDS =
     365          25 :         MEMDataset::Create("warp_temp", nXSize, nYSize, 0, GDT_Byte, nullptr);
     366             :     GDALRasterBandH hMEMBand =
     367          25 :         MEMCreateRasterBandEx(poMEMDS, 1, pabyPolyMask, GDT_Byte, 0, 0, false);
     368          25 :     poMEMDS->AddMEMBand(hMEMBand);
     369             : 
     370          25 :     GDALDatasetH hMemDS = GDALDataset::ToHandle(poMEMDS);
     371          25 :     double adfGeoTransform[6] = {0.0, 1.0, 0.0, 0.0, 0.0, 1.0};
     372          25 :     GDALSetGeoTransform(hMemDS, adfGeoTransform);
     373             : 
     374             :     /* -------------------------------------------------------------------- */
     375             :     /*      Burn the polygon into the mask with 1.0 values.                 */
     376             :     /* -------------------------------------------------------------------- */
     377          25 :     int nTargetBand = 1;
     378          25 :     double dfBurnValue = 255.0;
     379          25 :     char **papszRasterizeOptions = nullptr;
     380             : 
     381          25 :     if (CPLFetchBool(psWO->papszWarpOptions, "CUTLINE_ALL_TOUCHED", false))
     382             :         papszRasterizeOptions =
     383           5 :             CSLSetNameValue(papszRasterizeOptions, "ALL_TOUCHED", "TRUE");
     384             : 
     385          25 :     int anXYOff[2] = {nXOff, nYOff};
     386             : 
     387          25 :     CPLErr eErr = GDALRasterizeGeometries(
     388             :         hMemDS, 1, &nTargetBand, 1, &hPolygon, CutlineTransformer, anXYOff,
     389             :         &dfBurnValue, papszRasterizeOptions, nullptr, nullptr);
     390             : 
     391          25 :     CSLDestroy(papszRasterizeOptions);
     392             : 
     393             :     // Close and ensure data flushed to underlying array.
     394          25 :     GDALClose(hMemDS);
     395             : 
     396             :     /* -------------------------------------------------------------------- */
     397             :     /*      In the case with no blend distance, we just apply this as a     */
     398             :     /*      mask, zeroing out everything outside the polygon.               */
     399             :     /* -------------------------------------------------------------------- */
     400          25 :     if (psWO->dfCutlineBlendDist == 0.0)
     401             :     {
     402     2861310 :         for (int i = nXSize * nYSize - 1; i >= 0; i--)
     403             :         {
     404     2861290 :             if (pabyPolyMask[i] == 0)
     405     2602350 :                 static_cast<float *>(pValidityMask)[i] = 0.0;
     406             :         }
     407             :     }
     408             :     else
     409             :     {
     410           1 :         eErr = BlendMaskGenerator(nXOff, nYOff, nXSize, nYSize, pabyPolyMask,
     411             :                                   static_cast<float *>(pValidityMask), hPolygon,
     412             :                                   psWO->dfCutlineBlendDist);
     413             :     }
     414             : 
     415             :     /* -------------------------------------------------------------------- */
     416             :     /*      Clean up.                                                       */
     417             :     /* -------------------------------------------------------------------- */
     418          25 :     CPLFree(pabyPolyMask);
     419             : 
     420          25 :     return eErr;
     421             : }