LCOV - code coverage report
Current view: top level - alg - llrasterize.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 303 330 91.8 %
Date: 2025-12-05 02:43:06 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         326 : 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         326 :     if (!nPartCount)
      67             :     {
      68           0 :         return;
      69             :     }
      70             : 
      71         326 :     int n = 0;
      72         663 :     for (int part = 0; part < nPartCount; part++)
      73         337 :         n += panPartSize[part];
      74             : 
      75         652 :     std::vector<int> polyInts(n);
      76         652 :     std::vector<int> polyInts2;
      77         326 :     if (bAvoidBurningSamePoints)
      78          15 :         polyInts2.resize(n);
      79             : 
      80         326 :     double dminy = padfY[0];
      81         326 :     double dmaxy = padfY[0];
      82        9412 :     for (int i = 1; i < n; i++)
      83             :     {
      84        9086 :         if (padfY[i] < dminy)
      85             :         {
      86        1276 :             dminy = padfY[i];
      87             :         }
      88        7810 :         else if (padfY[i] > dmaxy)
      89             :         {
      90        2063 :             dmaxy = padfY[i];
      91             :         }
      92             :     }
      93         326 :     const int miny = static_cast<int>(std::max(0.0, dminy));
      94             :     const int maxy =
      95         326 :         static_cast<int>(std::min<double>(dmaxy, nRasterYSize - 1));
      96             : 
      97         326 :     constexpr int minx = 0;
      98         326 :     const int maxx = nRasterXSize - 1;
      99             : 
     100             :     // Fix in 1.3: count a vertex only once.
     101       12256 :     for (int y = miny; y <= maxy; y++)
     102             :     {
     103       11930 :         int partoffset = 0;
     104             : 
     105       11930 :         const double dy = y + 0.5;  // Center height of line.
     106             : 
     107       11930 :         int part = 0;
     108       11930 :         int ints = 0;
     109       11930 :         int ints2 = 0;
     110             : 
     111     7037230 :         for (int i = 0; i < n; i++)
     112             :         {
     113     7025300 :             if (i == partoffset + panPartSize[part])
     114             :             {
     115         150 :                 partoffset += panPartSize[part];
     116         150 :                 part++;
     117             :             }
     118             : 
     119     7025300 :             int ind1 = 0;
     120     7025300 :             int ind2 = 0;
     121     7025300 :             if (i == partoffset)
     122             :             {
     123       12080 :                 ind1 = partoffset + panPartSize[part] - 1;
     124       12080 :                 ind2 = partoffset;
     125             :             }
     126             :             else
     127             :             {
     128     7013220 :                 ind1 = i - 1;
     129     7013220 :                 ind2 = i;
     130             :             }
     131             : 
     132     7025300 :             double dy1 = padfY[ind1];
     133     7025300 :             double dy2 = padfY[ind2];
     134             : 
     135     7025300 :             if ((dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy))
     136     7001160 :                 continue;
     137             : 
     138       24262 :             double dx1 = 0.0;
     139       24262 :             double dx2 = 0.0;
     140       24262 :             if (dy1 < dy2)
     141             :             {
     142       12065 :                 dx1 = padfX[ind1];
     143       12065 :                 dx2 = padfX[ind2];
     144             :             }
     145       12197 :             else if (dy1 > dy2)
     146             :             {
     147       12075 :                 std::swap(dy1, dy2);
     148       12075 :                 dx2 = padfX[ind1];
     149       12075 :                 dx1 = padfX[ind2];
     150             :             }
     151             :             else  // if( fabs(dy1-dy2) < 1.e-6 )
     152             :             {
     153             : 
     154             :                 // AE: DO NOT skip bottom horizontal segments
     155             :                 // -Fill them separately-
     156         122 :                 if (padfX[ind1] > padfX[ind2])
     157             :                 {
     158          38 :                     const double dhorizontal_x1 = floor(padfX[ind2] + 0.5);
     159          38 :                     const double dhorizontal_x2 = floor(padfX[ind1] + 0.5);
     160             : 
     161          38 :                     if ((dhorizontal_x1 > static_cast<double>(maxx)) ||
     162          38 :                         (dhorizontal_x2 <= 0))
     163           0 :                         continue;
     164             :                     const int horizontal_x1 =
     165          38 :                         static_cast<int>(std::max(dhorizontal_x1, 0.0));
     166             :                     const int horizontal_x2 = static_cast<int>(
     167          38 :                         std::min<double>(dhorizontal_x2, nRasterXSize));
     168             : 
     169             :                     // Fill the horizontal segment (separately from the rest).
     170          38 :                     if (bAvoidBurningSamePoints)
     171             :                     {
     172          10 :                         polyInts2[ints2++] = horizontal_x1;
     173          10 :                         polyInts2[ints2++] = horizontal_x2;
     174             :                     }
     175             :                     else
     176             :                     {
     177          28 :                         pfnScanlineFunc(
     178             :                             pCBData, y, horizontal_x1, horizontal_x2 - 1,
     179             :                             dfVariant == nullptr ? 0 : dfVariant[0]);
     180             :                     }
     181             :                 }
     182             :                 // else: Skip top horizontal segments.
     183             :                 // They are already filled in the regular loop.
     184         122 :                 continue;
     185             :             }
     186             : 
     187       24140 :             if (dy < dy2 && dy >= dy1)
     188             :             {
     189             :                 const double intersect = std::clamp<double>(
     190       48092 :                     (dy - dy1) * (dx2 - dx1) / (dy2 - dy1) + dx1, INT_MIN,
     191       24046 :                     INT_MAX);
     192             : 
     193       24046 :                 polyInts[ints++] = static_cast<int>(floor(intersect + 0.5));
     194             :             }
     195             :         }
     196             : 
     197       11930 :         std::sort(polyInts.begin(), polyInts.begin() + ints);
     198       11930 :         std::sort(polyInts2.begin(), polyInts2.begin() + ints2);
     199             : 
     200       23953 :         for (int i = 0; i + 1 < ints; i += 2)
     201             :         {
     202       12023 :             if (polyInts[i] <= maxx && polyInts[i + 1] > minx)
     203             :             {
     204       11312 :                 pfnScanlineFunc(pCBData, y, polyInts[i], polyInts[i + 1] - 1,
     205             :                                 dfVariant == nullptr ? 0 : dfVariant[0]);
     206             :             }
     207             :         }
     208             : 
     209       11940 :         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 double dfX = padfX[i];
     241           5 :         const double dfY = padfY[i];
     242           5 :         double dfVariant = 0.0;
     243           5 :         if (padfVariant != nullptr)
     244           0 :             dfVariant = padfVariant[i];
     245             : 
     246           5 :         if (0 <= dfX && dfX < nRasterXSize && 0 <= dfY && dfY < nRasterYSize)
     247           5 :             pfnPointFunc(pCBData, static_cast<int>(dfY), static_cast<int>(dfX),
     248             :                          dfVariant);
     249             :     }
     250           5 : }
     251             : 
     252             : /************************************************************************/
     253             : /*                         GDALdllImageLine()                           */
     254             : /************************************************************************/
     255             : 
     256        1874 : void GDALdllImageLine(int nRasterXSize, int nRasterYSize, int nPartCount,
     257             :                       const int *panPartSize, const double *padfX,
     258             :                       const double *padfY, const double *padfVariant,
     259             :                       llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData)
     260             : {
     261        1874 :     if (!nPartCount)
     262           0 :         return;
     263             : 
     264        3748 :     for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
     265             :     {
     266       48667 :         for (int j = 1; j < panPartSize[i]; j++)
     267             :         {
     268       46793 :             double dfX = padfX[n + j - 1];
     269       46793 :             double dfY = padfY[n + j - 1];
     270       46793 :             double dfXEnd = padfX[n + j];
     271       46793 :             double dfYEnd = padfY[n + j];
     272             : 
     273             :             // Skip segments that are off the target region.
     274       46793 :             if ((dfY < 0.0 && dfYEnd < 0.0) ||
     275       46793 :                 (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
     276           0 :                 (dfX < 0.0 && dfXEnd < 0.0) ||
     277       46793 :                 (dfX > nRasterXSize && dfXEnd > nRasterXSize))
     278           0 :                 continue;
     279             : 
     280             :             // TODO: clamp coordinates to [0, nRasterXSize] * [0, nRasterYSize]
     281       46793 :             if (!(dfX >= INT_MIN && dfX <= INT_MAX && dfY >= INT_MIN &&
     282       46793 :                   dfY <= INT_MAX && dfXEnd >= INT_MIN && dfXEnd <= INT_MAX &&
     283       46793 :                   dfYEnd >= INT_MIN && dfYEnd <= INT_MAX))
     284             :             {
     285           0 :                 CPLErrorOnce(
     286             :                     CE_Warning, CPLE_AppDefined,
     287             :                     "GDALdllImageLine(): coordinates outside of int range");
     288           0 :                 continue;
     289             :             }
     290             : 
     291       46793 :             int iX = static_cast<int>(floor(dfX));
     292       46793 :             int iY = static_cast<int>(floor(dfY));
     293             : 
     294       46793 :             const int iX1 = static_cast<int>(floor(dfXEnd));
     295       46793 :             const int iY1 = static_cast<int>(floor(dfYEnd));
     296             : 
     297       46793 :             double dfVariant = 0.0;
     298       46793 :             double dfVariant1 = 0.0;
     299       46793 :             if (padfVariant != nullptr &&
     300       46758 :                 pCBData->eBurnValueSource != GBV_UserBurnValue)
     301             :             {
     302       46758 :                 dfVariant = padfVariant[n + j - 1];
     303       46758 :                 dfVariant1 = padfVariant[n + j];
     304             :             }
     305             : 
     306       46793 :             int nDeltaX = std::abs(iX1 - iX);
     307       46793 :             int nDeltaY = std::abs(iY1 - iY);
     308             : 
     309             :             // Step direction depends on line direction.
     310       46793 :             const int nXStep = (iX > iX1) ? -1 : 1;
     311       46793 :             const int nYStep = (iY > iY1) ? -1 : 1;
     312             : 
     313             :             // Determine the line slope.
     314       46793 :             if (nDeltaX >= nDeltaY)
     315             :             {
     316       38876 :                 const int nXError = nDeltaY << 1;
     317       38876 :                 const int nYError = nXError - (nDeltaX << 1);
     318       38876 :                 int nError = nXError - nDeltaX;
     319             :                 // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
     320             :                 // but if it is == 0, dfDeltaVariant is not really used, so any
     321             :                 // value is okay.
     322       38876 :                 const double dfDeltaVariant =
     323       38876 :                     nDeltaX == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaX;
     324             : 
     325             :                 // Do not burn the end point, unless we are in the last
     326             :                 // segment. This is to avoid burning twice intermediate points,
     327             :                 // which causes artifacts in Add mode
     328       38876 :                 if (j != panPartSize[i] - 1)
     329             :                 {
     330       37185 :                     nDeltaX--;
     331             :                 }
     332             : 
     333       67162 :                 while (nDeltaX-- >= 0)
     334             :                 {
     335       28286 :                     if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
     336             :                         iY < nRasterYSize)
     337       28089 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     338             : 
     339       28286 :                     dfVariant += dfDeltaVariant;
     340       28286 :                     iX += nXStep;
     341       28286 :                     if (nError > 0)
     342             :                     {
     343       11024 :                         iY += nYStep;
     344       11024 :                         nError += nYError;
     345             :                     }
     346             :                     else
     347             :                     {
     348       17262 :                         nError += nXError;
     349             :                     }
     350             :                 }
     351             :             }
     352             :             else
     353             :             {
     354        7917 :                 const int nXError = nDeltaX << 1;
     355        7917 :                 const int nYError = nXError - (nDeltaY << 1);
     356        7917 :                 int nError = nXError - nDeltaY;
     357             :                 // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
     358             :                 // but if it is == 0, dfDeltaVariant is not really used, so any
     359             :                 // value is okay.
     360        7917 :                 double dfDeltaVariant =
     361        7917 :                     nDeltaY == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaY;
     362             : 
     363             :                 // Do not burn the end point, unless we are in the last
     364             :                 // segment. This is to avoid burning twice intermediate points,
     365             :                 // which causes artifacts in Add mode
     366        7917 :                 if (j != panPartSize[i] - 1)
     367             :                 {
     368        7734 :                     nDeltaY--;
     369             :                 }
     370             : 
     371       16207 :                 while (nDeltaY-- >= 0)
     372             :                 {
     373        8290 :                     if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
     374             :                         iY < nRasterYSize)
     375        8201 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     376             : 
     377        8290 :                     dfVariant += dfDeltaVariant;
     378        8290 :                     iY += nYStep;
     379        8290 :                     if (nError > 0)
     380             :                     {
     381          35 :                         iX += nXStep;
     382          35 :                         nError += nYError;
     383             :                     }
     384             :                     else
     385             :                     {
     386        8255 :                         nError += nXError;
     387             :                     }
     388             :                 }
     389             :             }
     390             :         }
     391             :     }
     392             : }
     393             : 
     394             : /************************************************************************/
     395             : /*                     GDALdllImageLineAllTouched()                     */
     396             : /*                                                                      */
     397             : /*      This alternate line drawing algorithm attempts to ensure        */
     398             : /*      that every pixel touched at all by the line will get set.       */
     399             : /*      @param padfVariant should contain the values that are to be     */
     400             : /*      added to the burn value.  The values along the line between the */
     401             : /*      points will be linearly interpolated. These values are used     */
     402             : /*      only if pCBData->eBurnValueSource is set to something other     */
     403             : /*      than GBV_UserBurnValue. If NULL is passed, a monotonous line    */
     404             : /*      will be drawn with the burn value.                              */
     405             : /************************************************************************/
     406             : 
     407         141 : void GDALdllImageLineAllTouched(
     408             :     int nRasterXSize, int nRasterYSize, int nPartCount, const int *panPartSize,
     409             :     const double *padfX, const double *padfY, const double *padfVariant,
     410             :     llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData,
     411             :     bool bAvoidBurningSamePoints, bool bIntersectOnly)
     412             : 
     413             : {
     414             :     // This is an epsilon to detect geometries that are aligned with pixel
     415             :     // coordinates. Hard to find the right value. We put it to that value
     416             :     // to satisfy the scenarios of https://github.com/OSGeo/gdal/issues/7523
     417             :     // and https://github.com/OSGeo/gdal/issues/6414
     418         141 :     constexpr double EPSILON_INTERSECT_ONLY = 1e-4;
     419             : 
     420         141 :     if (!nPartCount)
     421           0 :         return;
     422             : 
     423         285 :     for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
     424             :     {
     425         288 :         std::set<std::pair<int, int>> lastBurntPoints;
     426         288 :         std::set<std::pair<int, int>> newBurntPoints;
     427             : 
     428        1000 :         for (int j = 1; j < panPartSize[i]; j++)
     429             :         {
     430         856 :             lastBurntPoints = std::move(newBurntPoints);
     431         856 :             newBurntPoints.clear();
     432             : 
     433         856 :             double dfX = padfX[n + j - 1];
     434         856 :             double dfY = padfY[n + j - 1];
     435             : 
     436         856 :             double dfXEnd = padfX[n + j];
     437         856 :             double dfYEnd = padfY[n + j];
     438             : 
     439             :             // Skip segments that are off the target region.
     440         856 :             if ((dfY < 0.0 && dfYEnd < 0.0) ||
     441         845 :                 (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
     442         831 :                 (dfX < 0.0 && dfXEnd < 0.0) ||
     443         800 :                 (dfX > nRasterXSize && dfXEnd > nRasterXSize))
     444         553 :                 continue;
     445             : 
     446             :             // TODO: clamp coordinates to [0, nRasterXSize] * [0, nRasterYSize]
     447         763 :             if (!(dfX >= INT_MIN && dfX <= INT_MAX && dfY >= INT_MIN &&
     448         763 :                   dfY <= INT_MAX && dfXEnd >= INT_MIN && dfXEnd <= INT_MAX &&
     449         763 :                   dfYEnd >= INT_MIN && dfYEnd <= INT_MAX))
     450             :             {
     451           0 :                 CPLErrorOnce(CE_Warning, CPLE_AppDefined,
     452             :                              "GDALdllImageLineAllTouched(): coordinates "
     453             :                              "outside of int range");
     454           0 :                 continue;
     455             :             }
     456             : 
     457         763 :             double dfVariant = 0.0;
     458         763 :             double dfVariantEnd = 0.0;
     459         763 :             if (padfVariant != nullptr &&
     460          11 :                 static_cast<GDALRasterizeInfo *>(pCBData)->eBurnValueSource !=
     461             :                     GBV_UserBurnValue)
     462             :             {
     463          11 :                 dfVariant = padfVariant[n + j - 1];
     464          11 :                 dfVariantEnd = padfVariant[n + j];
     465             :             }
     466             : 
     467             :             // Swap if needed so we can proceed from left2right (X increasing)
     468         763 :             if (dfX > dfXEnd)
     469             :             {
     470         272 :                 std::swap(dfX, dfXEnd);
     471         272 :                 std::swap(dfY, dfYEnd);
     472         272 :                 std::swap(dfVariant, dfVariantEnd);
     473             :             }
     474             : 
     475             :             // Special case for vertical lines.
     476             : 
     477         763 :             if (fabs(dfX - dfXEnd) < .01)
     478             :             {
     479         209 :                 if (bIntersectOnly)
     480             :                 {
     481         200 :                     if (std::abs(dfX - std::round(dfX)) <
     482         272 :                             EPSILON_INTERSECT_ONLY &&
     483          72 :                         std::abs(dfXEnd - std::round(dfXEnd)) <
     484             :                             EPSILON_INTERSECT_ONLY)
     485          71 :                         continue;
     486             :                 }
     487             : 
     488         138 :                 if (dfYEnd < dfY)
     489             :                 {
     490          73 :                     std::swap(dfY, dfYEnd);
     491          73 :                     std::swap(dfVariant, dfVariantEnd);
     492             :                 }
     493             : 
     494         138 :                 const int iX = static_cast<int>(floor(dfXEnd));
     495         138 :                 int iY = static_cast<int>(floor(dfY));
     496         138 :                 int iYEnd =
     497         138 :                     static_cast<int>(floor(dfYEnd - EPSILON_INTERSECT_ONLY));
     498             : 
     499         138 :                 if (iX < 0 || iX >= nRasterXSize)
     500           1 :                     continue;
     501             : 
     502         137 :                 double dfDeltaVariant = 0.0;
     503         137 :                 if (dfYEnd - dfY > 0.0)
     504         136 :                     dfDeltaVariant = (dfVariantEnd - dfVariant) /
     505         136 :                                      (dfYEnd - dfY);  // Per unit change in iY.
     506             : 
     507             :                 // Clip to the borders of the target region.
     508         137 :                 if (iY < 0)
     509           0 :                     iY = 0;
     510         137 :                 if (iYEnd >= nRasterYSize)
     511           0 :                     iYEnd = nRasterYSize - 1;
     512         137 :                 dfVariant += dfDeltaVariant * (iY - dfY);
     513             : 
     514         137 :                 if (padfVariant == nullptr)
     515             :                 {
     516        1289 :                     for (; iY <= iYEnd; iY++)
     517             :                     {
     518        1156 :                         if (bAvoidBurningSamePoints)
     519             :                         {
     520         146 :                             auto yx = std::pair<int, int>(iY, iX);
     521         146 :                             if (lastBurntPoints.find(yx) !=
     522         292 :                                 lastBurntPoints.end())
     523             :                             {
     524          13 :                                 continue;
     525             :                             }
     526         133 :                             newBurntPoints.insert(yx);
     527             :                         }
     528        1143 :                         pfnPointFunc(pCBData, iY, iX, 0.0);
     529             :                     }
     530             :                 }
     531             :                 else
     532             :                 {
     533          14 :                     for (; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant)
     534             :                     {
     535          10 :                         if (bAvoidBurningSamePoints)
     536             :                         {
     537           0 :                             auto yx = std::pair<int, int>(iY, iX);
     538           0 :                             if (lastBurntPoints.find(yx) !=
     539           0 :                                 lastBurntPoints.end())
     540             :                             {
     541           0 :                                 continue;
     542             :                             }
     543           0 :                             newBurntPoints.insert(yx);
     544             :                         }
     545          10 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     546             :                     }
     547             :                 }
     548             : 
     549         137 :                 continue;  // Next segment.
     550             :             }
     551             : 
     552         554 :             const double dfDeltaVariant =
     553         554 :                 (dfVariantEnd - dfVariant) /
     554         554 :                 (dfXEnd - dfX);  // Per unit change in iX.
     555             : 
     556             :             // Special case for horizontal lines.
     557         554 :             if (fabs(dfY - dfYEnd) < .01)
     558             :             {
     559         251 :                 if (bIntersectOnly)
     560             :                 {
     561         242 :                     if (std::abs(dfY - std::round(dfY)) <
     562         317 :                             EPSILON_INTERSECT_ONLY &&
     563          75 :                         std::abs(dfYEnd - std::round(dfYEnd)) <
     564             :                             EPSILON_INTERSECT_ONLY)
     565          74 :                         continue;
     566             :                 }
     567             : 
     568         177 :                 if (dfXEnd < dfX)
     569             :                 {
     570           0 :                     std::swap(dfX, dfXEnd);
     571           0 :                     std::swap(dfVariant, dfVariantEnd);
     572             :                 }
     573             : 
     574         177 :                 int iX = static_cast<int>(floor(dfX));
     575         177 :                 const int iY = static_cast<int>(floor(dfY));
     576         177 :                 int iXEnd =
     577         177 :                     static_cast<int>(floor(dfXEnd - EPSILON_INTERSECT_ONLY));
     578             : 
     579         177 :                 if (iY < 0 || iY >= nRasterYSize)
     580           0 :                     continue;
     581             : 
     582             :                 // Clip to the borders of the target region.
     583         177 :                 if (iX < 0)
     584           2 :                     iX = 0;
     585         177 :                 if (iXEnd >= nRasterXSize)
     586           1 :                     iXEnd = nRasterXSize - 1;
     587         177 :                 dfVariant += dfDeltaVariant * (iX - dfX);
     588             : 
     589         177 :                 if (padfVariant == nullptr)
     590             :                 {
     591        1163 :                     for (; iX <= iXEnd; iX++)
     592             :                     {
     593         990 :                         if (bAvoidBurningSamePoints)
     594             :                         {
     595         144 :                             auto yx = std::pair<int, int>(iY, iX);
     596         144 :                             if (lastBurntPoints.find(yx) !=
     597         288 :                                 lastBurntPoints.end())
     598             :                             {
     599          12 :                                 continue;
     600             :                             }
     601         132 :                             newBurntPoints.insert(yx);
     602             :                         }
     603         978 :                         pfnPointFunc(pCBData, iY, iX, 0.0);
     604             :                     }
     605             :                 }
     606             :                 else
     607             :                 {
     608          14 :                     for (; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant)
     609             :                     {
     610          10 :                         if (bAvoidBurningSamePoints)
     611             :                         {
     612           0 :                             auto yx = std::pair<int, int>(iY, iX);
     613           0 :                             if (lastBurntPoints.find(yx) !=
     614           0 :                                 lastBurntPoints.end())
     615             :                             {
     616           0 :                                 continue;
     617             :                             }
     618           0 :                             newBurntPoints.insert(yx);
     619             :                         }
     620          10 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     621             :                     }
     622             :                 }
     623             : 
     624         177 :                 continue;  // Next segment.
     625             :             }
     626             : 
     627             :             /* --------------------------------------------------------------------
     628             :              */
     629             :             /*      General case - left to right sloped. */
     630             :             /* --------------------------------------------------------------------
     631             :              */
     632         303 :             const double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
     633             : 
     634             :             // Clip segment in X.
     635         303 :             if (dfXEnd > nRasterXSize)
     636             :             {
     637          16 :                 dfYEnd -= (dfXEnd - nRasterXSize) * dfSlope;
     638          16 :                 dfXEnd = nRasterXSize;
     639             :             }
     640         303 :             if (dfX < 0.0)
     641             :             {
     642          18 :                 dfY += (0.0 - dfX) * dfSlope;
     643          18 :                 dfVariant += dfDeltaVariant * (0.0 - dfX);
     644          18 :                 dfX = 0.0;
     645             :             }
     646             : 
     647             :             // Clip segment in Y.
     648         303 :             if (dfYEnd > dfY)
     649             :             {
     650          74 :                 if (dfY < 0.0)
     651             :                 {
     652           7 :                     const double dfDiffX = (0.0 - dfY) / dfSlope;
     653           7 :                     dfX += dfDiffX;
     654           7 :                     dfVariant += dfDeltaVariant * dfDiffX;
     655           7 :                     dfY = 0.0;
     656             :                 }
     657          74 :                 if (dfYEnd >= nRasterYSize)
     658             :                 {
     659           9 :                     dfXEnd += (dfYEnd - nRasterYSize) / dfSlope;
     660           9 :                     if (dfXEnd > nRasterXSize)
     661           5 :                         dfXEnd = nRasterXSize;
     662             :                     // dfYEnd is no longer used afterwards, but for
     663             :                     // consistency it should be:
     664             :                     // dfYEnd = nRasterXSize;
     665             :                 }
     666             :             }
     667             :             else
     668             :             {
     669         229 :                 if (dfY >= nRasterYSize)
     670             :                 {
     671          13 :                     const double dfDiffX = (nRasterYSize - dfY) / dfSlope;
     672          13 :                     dfX += dfDiffX;
     673          13 :                     dfVariant += dfDeltaVariant * dfDiffX;
     674          13 :                     dfY = nRasterYSize;
     675             :                 }
     676         229 :                 if (dfYEnd < 0.0)
     677             :                 {
     678           7 :                     dfXEnd -= (dfYEnd - 0) / dfSlope;
     679             :                     // dfYEnd is no longer used afterwards, but for
     680             :                     // consistency it should be:
     681             :                     // dfYEnd = 0.0;
     682             :                 }
     683             :             }
     684             : 
     685             :             // Step from pixel to pixel.
     686        6738 :             while (dfX >= 0.0 && dfX < dfXEnd)
     687             :             {
     688        6435 :                 const int iX = static_cast<int>(floor(dfX));
     689        6435 :                 const int iY = static_cast<int>(floor(dfY));
     690             : 
     691             :                 // Burn in the current point.
     692             :                 // We should be able to drop the Y check because we clipped
     693             :                 // in Y, but there may be some error with all the small steps.
     694        6435 :                 if (iY >= 0 && iY < nRasterYSize)
     695             :                 {
     696        6336 :                     if (bAvoidBurningSamePoints)
     697             :                     {
     698          42 :                         auto yx = std::pair<int, int>(iY, iX);
     699          79 :                         if (lastBurntPoints.find(yx) == lastBurntPoints.end() &&
     700          79 :                             newBurntPoints.find(yx) == newBurntPoints.end())
     701             :                         {
     702          27 :                             newBurntPoints.insert(yx);
     703          27 :                             pfnPointFunc(pCBData, iY, iX, dfVariant);
     704             :                         }
     705             :                     }
     706             :                     else
     707             :                     {
     708        6294 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     709             :                     }
     710             :                 }
     711             : 
     712        6435 :                 double dfStepX = floor(dfX + 1.0) - dfX;
     713        6435 :                 double dfStepY = dfStepX * dfSlope;
     714             : 
     715             :                 // Step to right pixel without changing scanline?
     716        6435 :                 if (static_cast<int>(floor(dfY + dfStepY)) == iY)
     717             :                 {
     718        3840 :                     dfX += dfStepX;
     719        3840 :                     dfY += dfStepY;
     720        3840 :                     dfVariant += dfDeltaVariant * dfStepX;
     721             :                 }
     722        2595 :                 else if (dfSlope < 0)
     723             :                 {
     724        1780 :                     dfStepY = iY - dfY;
     725        1780 :                     if (dfStepY > -0.000000001)
     726         877 :                         dfStepY = -0.000000001;
     727             : 
     728        1780 :                     dfStepX = dfStepY / dfSlope;
     729        1780 :                     dfX += dfStepX;
     730        1780 :                     dfY += dfStepY;
     731        1780 :                     dfVariant += dfDeltaVariant * dfStepX;
     732             :                 }
     733             :                 else
     734             :                 {
     735         815 :                     dfStepY = (iY + 1) - dfY;
     736         815 :                     if (dfStepY < 0.000000001)
     737           1 :                         dfStepY = 0.000000001;
     738             : 
     739         815 :                     dfStepX = dfStepY / dfSlope;
     740         815 :                     dfX += dfStepX;
     741         815 :                     dfY += dfStepY;
     742         815 :                     dfVariant += dfDeltaVariant * dfStepX;
     743             :                 }
     744             :             }  // Next step along segment.
     745             :         }      // Next segment.
     746             :     }          // Next part.
     747             : }

Generated by: LCOV version 1.14