LCOV - code coverage report
Current view: top level - alg - llrasterize.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 269 314 85.7 %
Date: 2025-02-20 10:14:44 Functions: 4 4 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  GDAL
       4             :  * Purpose:  Vector polygon rasterization code.
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com>
       9             :  * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "gdal_alg_priv.h"
      16             : 
      17             : #include <cmath>
      18             : #include <cstdlib>
      19             : #include <cstring>
      20             : 
      21             : #include <algorithm>
      22             : #include <set>
      23             : #include <utility>
      24             : #include <vector>
      25             : 
      26             : #include "gdal_alg.h"
      27             : 
      28             : /************************************************************************/
      29             : /*                       dllImageFilledPolygon()                        */
      30             : /*                                                                      */
      31             : /*      Perform scanline conversion of the passed multi-ring            */
      32             : /*      polygon.  Note the polygon does not need to be explicitly       */
      33             : /*      closed.  The scanline function will be called with              */
      34             : /*      horizontal scanline chunks which may not be entirely            */
      35             : /*      contained within the valid raster area (in the X                */
      36             : /*      direction).                                                     */
      37             : /*                                                                      */
      38             : /*      NEW: Nodes' coordinate are kept as double  in order             */
      39             : /*      to compute accurately the intersections with the lines          */
      40             : /*                                                                      */
      41             : /*      A pixel is considered inside a polygon if its center            */
      42             : /*      falls inside the polygon. This is robust unless                 */
      43             : /*      the nodes are placed in the center of the pixels in which       */
      44             : /*      case, due to numerical inaccuracies, it's hard to predict       */
      45             : /*      if the pixel will be considered inside or outside the shape.    */
      46             : /************************************************************************/
      47             : 
      48             : /*
      49             :  * NOTE: This code was originally adapted from the gdImageFilledPolygon()
      50             :  * function in libgd.
      51             :  *
      52             :  * http://www.boutell.com/gd/
      53             :  *
      54             :  * It was later adapted for direct inclusion in GDAL and relicensed under
      55             :  * the GDAL MIT license (pulled from the OpenEV distribution).
      56             :  */
      57             : 
      58         105 : void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
      59             :                                int nPartCount, const int *panPartSize,
      60             :                                const double *padfX, const double *padfY,
      61             :                                const double *dfVariant,
      62             :                                llScanlineFunc pfnScanlineFunc,
      63             :                                GDALRasterizeInfo *pCBData,
      64             :                                bool bAvoidBurningSamePoints)
      65             : {
      66         105 :     if (!nPartCount)
      67             :     {
      68           0 :         return;
      69             :     }
      70             : 
      71         105 :     int n = 0;
      72         218 :     for (int part = 0; part < nPartCount; part++)
      73         113 :         n += panPartSize[part];
      74             : 
      75         210 :     std::vector<int> polyInts(n);
      76         210 :     std::vector<int> polyInts2;
      77         105 :     if (bAvoidBurningSamePoints)
      78          12 :         polyInts2.resize(n);
      79             : 
      80         105 :     double dminy = padfY[0];
      81         105 :     double dmaxy = padfY[0];
      82        6660 :     for (int i = 1; i < n; i++)
      83             :     {
      84        6555 :         if (padfY[i] < dminy)
      85             :         {
      86         907 :             dminy = padfY[i];
      87             :         }
      88        6555 :         if (padfY[i] > dmaxy)
      89             :         {
      90        1322 :             dmaxy = padfY[i];
      91             :         }
      92             :     }
      93         105 :     int miny = static_cast<int>(dminy);
      94         105 :     int maxy = static_cast<int>(dmaxy);
      95             : 
      96         105 :     if (miny < 0)
      97           9 :         miny = 0;
      98         105 :     if (maxy >= nRasterYSize)
      99          15 :         maxy = nRasterYSize - 1;
     100             : 
     101         105 :     int minx = 0;
     102         105 :     const int maxx = nRasterXSize - 1;
     103             : 
     104             :     // Fix in 1.3: count a vertex only once.
     105        9248 :     for (int y = miny; y <= maxy; y++)
     106             :     {
     107        9143 :         int partoffset = 0;
     108             : 
     109        9143 :         const double dy = y + 0.5;  // Center height of line.
     110             : 
     111        9143 :         int part = 0;
     112        9143 :         int ints = 0;
     113        9143 :         int ints2 = 0;
     114             : 
     115     6964140 :         for (int i = 0; i < n; i++)
     116             :         {
     117     6955000 :             if (i == partoffset + panPartSize[part])
     118             :             {
     119         114 :                 partoffset += panPartSize[part];
     120         114 :                 part++;
     121             :             }
     122             : 
     123     6955000 :             int ind1 = 0;
     124     6955000 :             int ind2 = 0;
     125     6955000 :             if (i == partoffset)
     126             :             {
     127        9257 :                 ind1 = partoffset + panPartSize[part] - 1;
     128        9257 :                 ind2 = partoffset;
     129             :             }
     130             :             else
     131             :             {
     132     6945740 :                 ind1 = i - 1;
     133     6945740 :                 ind2 = i;
     134             :             }
     135             : 
     136     6955000 :             double dy1 = padfY[ind1];
     137     6955000 :             double dy2 = padfY[ind2];
     138             : 
     139     6955000 :             if ((dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy))
     140     6936690 :                 continue;
     141             : 
     142       18394 :             double dx1 = 0.0;
     143       18394 :             double dx2 = 0.0;
     144       18394 :             if (dy1 < dy2)
     145             :             {
     146        9157 :                 dx1 = padfX[ind1];
     147        9157 :                 dx2 = padfX[ind2];
     148             :             }
     149        9237 :             else if (dy1 > dy2)
     150             :             {
     151        9157 :                 std::swap(dy1, dy2);
     152        9157 :                 dx2 = padfX[ind1];
     153        9157 :                 dx1 = padfX[ind2];
     154             :             }
     155             :             else  // if( fabs(dy1-dy2) < 1.e-6 )
     156             :             {
     157             : 
     158             :                 // AE: DO NOT skip bottom horizontal segments
     159             :                 // -Fill them separately-
     160          80 :                 if (padfX[ind1] > padfX[ind2])
     161             :                 {
     162          24 :                     const int horizontal_x1 =
     163          24 :                         static_cast<int>(floor(padfX[ind2] + 0.5));
     164          24 :                     const int horizontal_x2 =
     165          24 :                         static_cast<int>(floor(padfX[ind1] + 0.5));
     166             : 
     167          24 :                     if ((horizontal_x1 > maxx) || (horizontal_x2 <= minx))
     168           0 :                         continue;
     169             : 
     170             :                     // Fill the horizontal segment (separately from the rest).
     171          24 :                     if (bAvoidBurningSamePoints)
     172             :                     {
     173          10 :                         polyInts2[ints2++] = horizontal_x1;
     174          10 :                         polyInts2[ints2++] = horizontal_x2;
     175             :                     }
     176             :                     else
     177             :                     {
     178          14 :                         pfnScanlineFunc(
     179             :                             pCBData, y, horizontal_x1, horizontal_x2 - 1,
     180             :                             dfVariant == nullptr ? 0 : dfVariant[0]);
     181             :                     }
     182             :                 }
     183             :                 // else: Skip top horizontal segments.
     184             :                 // They are already filled in the regular loop.
     185          80 :                 continue;
     186             :             }
     187             : 
     188       18314 :             if (dy < dy2 && dy >= dy1)
     189             :             {
     190       18258 :                 const double intersect =
     191       18258 :                     (dy - dy1) * (dx2 - dx1) / (dy2 - dy1) + dx1;
     192             : 
     193       18258 :                 polyInts[ints++] = static_cast<int>(floor(intersect + 0.5));
     194             :             }
     195             :         }
     196             : 
     197        9143 :         std::sort(polyInts.begin(), polyInts.begin() + ints);
     198        9143 :         std::sort(polyInts2.begin(), polyInts2.begin() + ints2);
     199             : 
     200       18272 :         for (int i = 0; i + 1 < ints; i += 2)
     201             :         {
     202        9129 :             if (polyInts[i] <= maxx && polyInts[i + 1] > minx)
     203             :             {
     204        8461 :                 pfnScanlineFunc(pCBData, y, polyInts[i], polyInts[i + 1] - 1,
     205             :                                 dfVariant == nullptr ? 0 : dfVariant[0]);
     206             :             }
     207             :         }
     208             : 
     209        9153 :         for (int i2 = 0, i = 0; i2 + 1 < ints2; i2 += 2)
     210             :         {
     211          10 :             if (polyInts2[i2] <= maxx && polyInts2[i2 + 1] > minx)
     212             :             {
     213             :                 // "synchronize" polyInts[i] with polyInts2[i2]
     214          10 :                 while (i + 1 < ints && polyInts[i] < polyInts2[i2])
     215           0 :                     i += 2;
     216             :                 // Only burn if we don't have a common segment between
     217             :                 // polyInts[] and polyInts2[]
     218          10 :                 if (i + 1 >= ints || polyInts[i] != polyInts2[i2])
     219             :                 {
     220           4 :                     pfnScanlineFunc(pCBData, y, polyInts2[i2],
     221           2 :                                     polyInts2[i2 + 1] - 1,
     222             :                                     dfVariant == nullptr ? 0 : dfVariant[0]);
     223             :                 }
     224             :             }
     225             :         }
     226             :     }
     227             : }
     228             : 
     229             : /************************************************************************/
     230             : /*                         GDALdllImagePoint()                          */
     231             : /************************************************************************/
     232             : 
     233           5 : void GDALdllImagePoint(int nRasterXSize, int nRasterYSize, int nPartCount,
     234             :                        const int * /*panPartSize*/, const double *padfX,
     235             :                        const double *padfY, const double *padfVariant,
     236             :                        llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData)
     237             : {
     238          10 :     for (int i = 0; i < nPartCount; i++)
     239             :     {
     240           5 :         const int nX = static_cast<int>(floor(padfX[i]));
     241           5 :         const int nY = static_cast<int>(floor(padfY[i]));
     242           5 :         double dfVariant = 0.0;
     243           5 :         if (padfVariant != nullptr)
     244           0 :             dfVariant = padfVariant[i];
     245             : 
     246           5 :         if (0 <= nX && nX < nRasterXSize && 0 <= nY && nY < nRasterYSize)
     247           5 :             pfnPointFunc(pCBData, nY, nX, dfVariant);
     248             :     }
     249           5 : }
     250             : 
     251             : /************************************************************************/
     252             : /*                         GDALdllImageLine()                           */
     253             : /************************************************************************/
     254             : 
     255        1874 : void GDALdllImageLine(int nRasterXSize, int nRasterYSize, int nPartCount,
     256             :                       const int *panPartSize, const double *padfX,
     257             :                       const double *padfY, const double *padfVariant,
     258             :                       llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData)
     259             : {
     260        1874 :     if (!nPartCount)
     261           0 :         return;
     262             : 
     263        3748 :     for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
     264             :     {
     265       48667 :         for (int j = 1; j < panPartSize[i]; j++)
     266             :         {
     267       46793 :             int iX = static_cast<int>(floor(padfX[n + j - 1]));
     268       46793 :             int iY = static_cast<int>(floor(padfY[n + j - 1]));
     269             : 
     270       46793 :             const int iX1 = static_cast<int>(floor(padfX[n + j]));
     271       46793 :             const int iY1 = static_cast<int>(floor(padfY[n + j]));
     272             : 
     273       46793 :             double dfVariant = 0.0;
     274       46793 :             double dfVariant1 = 0.0;
     275       46793 :             if (padfVariant != nullptr &&
     276       46758 :                 pCBData->eBurnValueSource != GBV_UserBurnValue)
     277             :             {
     278       46758 :                 dfVariant = padfVariant[n + j - 1];
     279       46758 :                 dfVariant1 = padfVariant[n + j];
     280             :             }
     281             : 
     282       46793 :             int nDeltaX = std::abs(iX1 - iX);
     283       46793 :             int nDeltaY = std::abs(iY1 - iY);
     284             : 
     285             :             // Step direction depends on line direction.
     286       46793 :             const int nXStep = (iX > iX1) ? -1 : 1;
     287       46793 :             const int nYStep = (iY > iY1) ? -1 : 1;
     288             : 
     289             :             // Determine the line slope.
     290       46793 :             if (nDeltaX >= nDeltaY)
     291             :             {
     292       38876 :                 const int nXError = nDeltaY << 1;
     293       38876 :                 const int nYError = nXError - (nDeltaX << 1);
     294       38876 :                 int nError = nXError - nDeltaX;
     295             :                 // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
     296             :                 // but if it is == 0, dfDeltaVariant is not really used, so any
     297             :                 // value is okay.
     298       38876 :                 const double dfDeltaVariant =
     299       38876 :                     nDeltaX == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaX;
     300             : 
     301             :                 // Do not burn the end point, unless we are in the last
     302             :                 // segment. This is to avoid burning twice intermediate points,
     303             :                 // which causes artifacts in Add mode
     304       38876 :                 if (j != panPartSize[i] - 1)
     305             :                 {
     306       37185 :                     nDeltaX--;
     307             :                 }
     308             : 
     309       67162 :                 while (nDeltaX-- >= 0)
     310             :                 {
     311       28286 :                     if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
     312             :                         iY < nRasterYSize)
     313       28089 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     314             : 
     315       28286 :                     dfVariant += dfDeltaVariant;
     316       28286 :                     iX += nXStep;
     317       28286 :                     if (nError > 0)
     318             :                     {
     319       11024 :                         iY += nYStep;
     320       11024 :                         nError += nYError;
     321             :                     }
     322             :                     else
     323             :                     {
     324       17262 :                         nError += nXError;
     325             :                     }
     326             :                 }
     327             :             }
     328             :             else
     329             :             {
     330        7917 :                 const int nXError = nDeltaX << 1;
     331        7917 :                 const int nYError = nXError - (nDeltaY << 1);
     332        7917 :                 int nError = nXError - nDeltaY;
     333             :                 // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
     334             :                 // but if it is == 0, dfDeltaVariant is not really used, so any
     335             :                 // value is okay.
     336        7917 :                 double dfDeltaVariant =
     337        7917 :                     nDeltaY == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaY;
     338             : 
     339             :                 // Do not burn the end point, unless we are in the last
     340             :                 // segment. This is to avoid burning twice intermediate points,
     341             :                 // which causes artifacts in Add mode
     342        7917 :                 if (j != panPartSize[i] - 1)
     343             :                 {
     344        7734 :                     nDeltaY--;
     345             :                 }
     346             : 
     347       16207 :                 while (nDeltaY-- >= 0)
     348             :                 {
     349        8290 :                     if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
     350             :                         iY < nRasterYSize)
     351        8201 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     352             : 
     353        8290 :                     dfVariant += dfDeltaVariant;
     354        8290 :                     iY += nYStep;
     355        8290 :                     if (nError > 0)
     356             :                     {
     357          35 :                         iX += nXStep;
     358          35 :                         nError += nYError;
     359             :                     }
     360             :                     else
     361             :                     {
     362        8255 :                         nError += nXError;
     363             :                     }
     364             :                 }
     365             :             }
     366             :         }
     367             :     }
     368             : }
     369             : 
     370             : /************************************************************************/
     371             : /*                     GDALdllImageLineAllTouched()                     */
     372             : /*                                                                      */
     373             : /*      This alternate line drawing algorithm attempts to ensure        */
     374             : /*      that every pixel touched at all by the line will get set.       */
     375             : /*      @param padfVariant should contain the values that are to be     */
     376             : /*      added to the burn value.  The values along the line between the */
     377             : /*      points will be linearly interpolated. These values are used     */
     378             : /*      only if pCBData->eBurnValueSource is set to something other     */
     379             : /*      than GBV_UserBurnValue. If NULL is passed, a monotonous line    */
     380             : /*      will be drawn with the burn value.                              */
     381             : /************************************************************************/
     382             : 
     383          39 : void GDALdllImageLineAllTouched(
     384             :     int nRasterXSize, int nRasterYSize, int nPartCount, const int *panPartSize,
     385             :     const double *padfX, const double *padfY, const double *padfVariant,
     386             :     llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData,
     387             :     bool bAvoidBurningSamePoints, bool bIntersectOnly)
     388             : 
     389             : {
     390             :     // This is an epsilon to detect geometries that are aligned with pixel
     391             :     // coordinates. Hard to find the right value. We put it to that value
     392             :     // to satisfy the scenarios of https://github.com/OSGeo/gdal/issues/7523
     393             :     // and https://github.com/OSGeo/gdal/issues/6414
     394          39 :     constexpr double EPSILON_INTERSECT_ONLY = 1e-4;
     395             : 
     396          39 :     if (!nPartCount)
     397           0 :         return;
     398             : 
     399          78 :     for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
     400             :     {
     401          78 :         std::set<std::pair<int, int>> lastBurntPoints;
     402          78 :         std::set<std::pair<int, int>> newBurntPoints;
     403             : 
     404         186 :         for (int j = 1; j < panPartSize[i]; j++)
     405             :         {
     406         147 :             lastBurntPoints = std::move(newBurntPoints);
     407         147 :             newBurntPoints.clear();
     408             : 
     409         147 :             double dfX = padfX[n + j - 1];
     410         147 :             double dfY = padfY[n + j - 1];
     411             : 
     412         147 :             double dfXEnd = padfX[n + j];
     413         147 :             double dfYEnd = padfY[n + j];
     414             : 
     415         147 :             double dfVariant = 0.0;
     416         147 :             double dfVariantEnd = 0.0;
     417         147 :             if (padfVariant != nullptr &&
     418           0 :                 static_cast<GDALRasterizeInfo *>(pCBData)->eBurnValueSource !=
     419             :                     GBV_UserBurnValue)
     420             :             {
     421           0 :                 dfVariant = padfVariant[n + j - 1];
     422           0 :                 dfVariantEnd = padfVariant[n + j];
     423             :             }
     424             : 
     425             :             // Skip segments that are off the target region.
     426         147 :             if ((dfY < 0.0 && dfYEnd < 0.0) ||
     427         146 :                 (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
     428         145 :                 (dfX < 0.0 && dfXEnd < 0.0) ||
     429         144 :                 (dfX > nRasterXSize && dfXEnd > nRasterXSize))
     430         130 :                 continue;
     431             : 
     432             :             // Swap if needed so we can proceed from left2right (X increasing)
     433         143 :             if (dfX > dfXEnd)
     434             :             {
     435          43 :                 std::swap(dfX, dfXEnd);
     436          43 :                 std::swap(dfY, dfYEnd);
     437          43 :                 std::swap(dfVariant, dfVariantEnd);
     438             :             }
     439             : 
     440             :             // Special case for vertical lines.
     441         143 :             if (floor(dfX) == floor(dfXEnd) || fabs(dfX - dfXEnd) < .01)
     442             :             {
     443          65 :                 if (bIntersectOnly)
     444             :                 {
     445          56 :                     if (std::abs(dfX - std::round(dfX)) <
     446          67 :                             EPSILON_INTERSECT_ONLY &&
     447          11 :                         std::abs(dfXEnd - std::round(dfXEnd)) <
     448             :                             EPSILON_INTERSECT_ONLY)
     449          11 :                         continue;
     450             :                 }
     451             : 
     452          54 :                 if (dfYEnd < dfY)
     453             :                 {
     454          29 :                     std::swap(dfY, dfYEnd);
     455          29 :                     std::swap(dfVariant, dfVariantEnd);
     456             :                 }
     457             : 
     458          54 :                 const int iX = static_cast<int>(floor(dfXEnd));
     459          54 :                 int iY = static_cast<int>(floor(dfY));
     460          54 :                 int iYEnd =
     461          54 :                     static_cast<int>(floor(dfYEnd - EPSILON_INTERSECT_ONLY));
     462             : 
     463          54 :                 if (iX < 0 || iX >= nRasterXSize)
     464           0 :                     continue;
     465             : 
     466          54 :                 double dfDeltaVariant = 0.0;
     467          54 :                 if (dfYEnd - dfY > 0.0)
     468          54 :                     dfDeltaVariant = (dfVariantEnd - dfVariant) /
     469          54 :                                      (dfYEnd - dfY);  // Per unit change in iY.
     470             : 
     471             :                 // Clip to the borders of the target region.
     472          54 :                 if (iY < 0)
     473           0 :                     iY = 0;
     474          54 :                 if (iYEnd >= nRasterYSize)
     475           0 :                     iYEnd = nRasterYSize - 1;
     476          54 :                 dfVariant += dfDeltaVariant * (iY - dfY);
     477             : 
     478          54 :                 if (padfVariant == nullptr)
     479             :                 {
     480         953 :                     for (; iY <= iYEnd; iY++)
     481             :                     {
     482         899 :                         if (bAvoidBurningSamePoints)
     483             :                         {
     484         136 :                             auto yx = std::pair<int, int>(iY, iX);
     485         136 :                             if (lastBurntPoints.find(yx) !=
     486         272 :                                 lastBurntPoints.end())
     487             :                             {
     488          11 :                                 continue;
     489             :                             }
     490         125 :                             newBurntPoints.insert(yx);
     491             :                         }
     492         888 :                         pfnPointFunc(pCBData, iY, iX, 0.0);
     493             :                     }
     494             :                 }
     495             :                 else
     496             :                 {
     497           0 :                     for (; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant)
     498             :                     {
     499           0 :                         if (bAvoidBurningSamePoints)
     500             :                         {
     501           0 :                             auto yx = std::pair<int, int>(iY, iX);
     502           0 :                             if (lastBurntPoints.find(yx) !=
     503           0 :                                 lastBurntPoints.end())
     504             :                             {
     505           0 :                                 continue;
     506             :                             }
     507           0 :                             newBurntPoints.insert(yx);
     508             :                         }
     509           0 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     510             :                     }
     511             :                 }
     512             : 
     513          54 :                 continue;  // Next segment.
     514             :             }
     515             : 
     516          78 :             const double dfDeltaVariant =
     517          78 :                 (dfVariantEnd - dfVariant) /
     518          78 :                 (dfXEnd - dfX);  // Per unit change in iX.
     519             : 
     520             :             // Special case for horizontal lines.
     521          78 :             if (floor(dfY) == floor(dfYEnd) || fabs(dfY - dfYEnd) < .01)
     522             :             {
     523          61 :                 if (bIntersectOnly)
     524             :                 {
     525          52 :                     if (std::abs(dfY - std::round(dfY)) <
     526          59 :                             EPSILON_INTERSECT_ONLY &&
     527           7 :                         std::abs(dfYEnd - std::round(dfYEnd)) <
     528             :                             EPSILON_INTERSECT_ONLY)
     529           7 :                         continue;
     530             :                 }
     531             : 
     532          54 :                 if (dfXEnd < dfX)
     533             :                 {
     534           0 :                     std::swap(dfX, dfXEnd);
     535           0 :                     std::swap(dfVariant, dfVariantEnd);
     536             :                 }
     537             : 
     538          54 :                 int iX = static_cast<int>(floor(dfX));
     539          54 :                 const int iY = static_cast<int>(floor(dfY));
     540          54 :                 int iXEnd =
     541          54 :                     static_cast<int>(floor(dfXEnd - EPSILON_INTERSECT_ONLY));
     542             : 
     543          54 :                 if (iY < 0 || iY >= nRasterYSize)
     544           0 :                     continue;
     545             : 
     546             :                 // Clip to the borders of the target region.
     547          54 :                 if (iX < 0)
     548           0 :                     iX = 0;
     549          54 :                 if (iXEnd >= nRasterXSize)
     550           0 :                     iXEnd = nRasterXSize - 1;
     551          54 :                 dfVariant += dfDeltaVariant * (iX - dfX);
     552             : 
     553          54 :                 if (padfVariant == nullptr)
     554             :                 {
     555         745 :                     for (; iX <= iXEnd; iX++)
     556             :                     {
     557         691 :                         if (bAvoidBurningSamePoints)
     558             :                         {
     559         134 :                             auto yx = std::pair<int, int>(iY, iX);
     560         134 :                             if (lastBurntPoints.find(yx) !=
     561         268 :                                 lastBurntPoints.end())
     562             :                             {
     563           8 :                                 continue;
     564             :                             }
     565         126 :                             newBurntPoints.insert(yx);
     566             :                         }
     567         683 :                         pfnPointFunc(pCBData, iY, iX, 0.0);
     568             :                     }
     569             :                 }
     570             :                 else
     571             :                 {
     572           0 :                     for (; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant)
     573             :                     {
     574           0 :                         if (bAvoidBurningSamePoints)
     575             :                         {
     576           0 :                             auto yx = std::pair<int, int>(iY, iX);
     577           0 :                             if (lastBurntPoints.find(yx) !=
     578           0 :                                 lastBurntPoints.end())
     579             :                             {
     580           0 :                                 continue;
     581             :                             }
     582           0 :                             newBurntPoints.insert(yx);
     583             :                         }
     584           0 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     585             :                     }
     586             :                 }
     587             : 
     588          54 :                 continue;  // Next segment.
     589             :             }
     590             : 
     591             :             /* --------------------------------------------------------------------
     592             :              */
     593             :             /*      General case - left to right sloped. */
     594             :             /* --------------------------------------------------------------------
     595             :              */
     596          17 :             const double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
     597             : 
     598             :             // Clip segment in X.
     599          17 :             if (dfXEnd > nRasterXSize)
     600             :             {
     601           0 :                 dfYEnd -= (dfXEnd - nRasterXSize) * dfSlope;
     602           0 :                 dfXEnd = nRasterXSize;
     603             :             }
     604          17 :             if (dfX < 0.0)
     605             :             {
     606           1 :                 dfY += (0.0 - dfX) * dfSlope;
     607           1 :                 dfVariant += dfDeltaVariant * (0.0 - dfX);
     608           1 :                 dfX = 0.0;
     609             :             }
     610             : 
     611             :             // Clip segment in Y.
     612          17 :             if (dfYEnd > dfY)
     613             :             {
     614           6 :                 if (dfY < 0.0)
     615             :                 {
     616           0 :                     const double dfDiffX = (0.0 - dfY) / dfSlope;
     617           0 :                     dfX += dfDiffX;
     618           0 :                     dfVariant += dfDeltaVariant * dfDiffX;
     619           0 :                     dfY = 0.0;
     620             :                 }
     621           6 :                 if (dfYEnd >= nRasterYSize)
     622             :                 {
     623           2 :                     dfXEnd += (dfYEnd - nRasterYSize) / dfSlope;
     624           2 :                     if (dfXEnd > nRasterXSize)
     625           2 :                         dfXEnd = nRasterXSize;
     626             :                     // dfYEnd is no longer used afterwards, but for
     627             :                     // consistency it should be:
     628             :                     // dfYEnd = nRasterXSize;
     629             :                 }
     630             :             }
     631             :             else
     632             :             {
     633          11 :                 if (dfY >= nRasterYSize)
     634             :                 {
     635           0 :                     const double dfDiffX = (nRasterYSize - dfY) / dfSlope;
     636           0 :                     dfX += dfDiffX;
     637           0 :                     dfVariant += dfDeltaVariant * dfDiffX;
     638           0 :                     dfY = nRasterYSize;
     639             :                 }
     640          11 :                 if (dfYEnd < 0.0)
     641             :                 {
     642           0 :                     dfXEnd -= (dfYEnd - 0) / dfSlope;
     643             :                     // dfYEnd is no longer used afterwards, but for
     644             :                     // consistency it should be:
     645             :                     // dfYEnd = 0.0;
     646             :                 }
     647             :             }
     648             : 
     649             :             // Step from pixel to pixel.
     650        4411 :             while (dfX >= 0.0 && dfX < dfXEnd)
     651             :             {
     652        4394 :                 const int iX = static_cast<int>(floor(dfX));
     653        4394 :                 const int iY = static_cast<int>(floor(dfY));
     654             : 
     655             :                 // Burn in the current point.
     656             :                 // We should be able to drop the Y check because we clipped
     657             :                 // in Y, but there may be some error with all the small steps.
     658        4394 :                 if (iY >= 0 && iY < nRasterYSize)
     659             :                 {
     660        4393 :                     if (bAvoidBurningSamePoints)
     661             :                     {
     662          29 :                         auto yx = std::pair<int, int>(iY, iX);
     663          57 :                         if (lastBurntPoints.find(yx) == lastBurntPoints.end() &&
     664          57 :                             newBurntPoints.find(yx) == newBurntPoints.end())
     665             :                         {
     666          19 :                             newBurntPoints.insert(yx);
     667          19 :                             pfnPointFunc(pCBData, iY, iX, dfVariant);
     668             :                         }
     669             :                     }
     670             :                     else
     671             :                     {
     672        4364 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     673             :                     }
     674             :                 }
     675             : 
     676        4394 :                 double dfStepX = floor(dfX + 1.0) - dfX;
     677        4394 :                 double dfStepY = dfStepX * dfSlope;
     678             : 
     679             :                 // Step to right pixel without changing scanline?
     680        4394 :                 if (static_cast<int>(floor(dfY + dfStepY)) == iY)
     681             :                 {
     682        2907 :                     dfX += dfStepX;
     683        2907 :                     dfY += dfStepY;
     684        2907 :                     dfVariant += dfDeltaVariant * dfStepX;
     685             :                 }
     686        1487 :                 else if (dfSlope < 0)
     687             :                 {
     688         956 :                     dfStepY = iY - dfY;
     689         956 :                     if (dfStepY > -0.000000001)
     690         478 :                         dfStepY = -0.000000001;
     691             : 
     692         956 :                     dfStepX = dfStepY / dfSlope;
     693         956 :                     dfX += dfStepX;
     694         956 :                     dfY += dfStepY;
     695         956 :                     dfVariant += dfDeltaVariant * dfStepX;
     696             :                 }
     697             :                 else
     698             :                 {
     699         531 :                     dfStepY = (iY + 1) - dfY;
     700         531 :                     if (dfStepY < 0.000000001)
     701           0 :                         dfStepY = 0.000000001;
     702             : 
     703         531 :                     dfStepX = dfStepY / dfSlope;
     704         531 :                     dfX += dfStepX;
     705         531 :                     dfY += dfStepY;
     706         531 :                     dfVariant += dfDeltaVariant * dfStepX;
     707             :                 }
     708             :             }  // Next step along segment.
     709             :         }      // Next segment.
     710             :     }          // Next part.
     711             : }

Generated by: LCOV version 1.14