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