LCOV - code coverage report
Current view: top level - alg - llrasterize.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 286 314 91.1 %
Date: 2025-05-20 14:45:53 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         159 : 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         159 :     if (!nPartCount)
      67             :     {
      68           0 :         return;
      69             :     }
      70             : 
      71         159 :     int n = 0;
      72         329 :     for (int part = 0; part < nPartCount; part++)
      73         170 :         n += panPartSize[part];
      74             : 
      75         318 :     std::vector<int> polyInts(n);
      76         318 :     std::vector<int> polyInts2;
      77         159 :     if (bAvoidBurningSamePoints)
      78          15 :         polyInts2.resize(n);
      79             : 
      80         159 :     double dminy = padfY[0];
      81         159 :     double dmaxy = padfY[0];
      82        7203 :     for (int i = 1; i < n; i++)
      83             :     {
      84        7044 :         if (padfY[i] < dminy)
      85             :         {
      86        1015 :             dminy = padfY[i];
      87             :         }
      88        7044 :         if (padfY[i] > dmaxy)
      89             :         {
      90        1362 :             dmaxy = padfY[i];
      91             :         }
      92             :     }
      93         159 :     int miny = static_cast<int>(dminy);
      94         159 :     int maxy = static_cast<int>(dmaxy);
      95             : 
      96         159 :     if (miny < 0)
      97          10 :         miny = 0;
      98         159 :     if (maxy >= nRasterYSize)
      99          19 :         maxy = nRasterYSize - 1;
     100             : 
     101         159 :     int minx = 0;
     102         159 :     const int maxx = nRasterXSize - 1;
     103             : 
     104             :     // Fix in 1.3: count a vertex only once.
     105        9584 :     for (int y = miny; y <= maxy; y++)
     106             :     {
     107        9425 :         int partoffset = 0;
     108             : 
     109        9425 :         const double dy = y + 0.5;  // Center height of line.
     110             : 
     111        9425 :         int part = 0;
     112        9425 :         int ints = 0;
     113        9425 :         int ints2 = 0;
     114             : 
     115     6971380 :         for (int i = 0; i < n; i++)
     116             :         {
     117     6961950 :             if (i == partoffset + panPartSize[part])
     118             :             {
     119         150 :                 partoffset += panPartSize[part];
     120         150 :                 part++;
     121             :             }
     122             : 
     123     6961950 :             int ind1 = 0;
     124     6961950 :             int ind2 = 0;
     125     6961950 :             if (i == partoffset)
     126             :             {
     127        9575 :                 ind1 = partoffset + panPartSize[part] - 1;
     128        9575 :                 ind2 = partoffset;
     129             :             }
     130             :             else
     131             :             {
     132     6952380 :                 ind1 = i - 1;
     133     6952380 :                 ind2 = i;
     134             :             }
     135             : 
     136     6961950 :             double dy1 = padfY[ind1];
     137     6961950 :             double dy2 = padfY[ind2];
     138             : 
     139     6961950 :             if ((dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy))
     140     6943100 :                 continue;
     141             : 
     142       18962 :             double dx1 = 0.0;
     143       18962 :             double dx2 = 0.0;
     144       18962 :             if (dy1 < dy2)
     145             :             {
     146        9423 :                 dx1 = padfX[ind1];
     147        9423 :                 dx2 = padfX[ind2];
     148             :             }
     149        9539 :             else if (dy1 > dy2)
     150             :             {
     151        9423 :                 std::swap(dy1, dy2);
     152        9423 :                 dx2 = padfX[ind1];
     153        9423 :                 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         116 :                 if (padfX[ind1] > padfX[ind2])
     161             :                 {
     162          36 :                     const int horizontal_x1 =
     163          36 :                         static_cast<int>(floor(padfX[ind2] + 0.5));
     164          36 :                     const int horizontal_x2 =
     165          36 :                         static_cast<int>(floor(padfX[ind1] + 0.5));
     166             : 
     167          36 :                     if ((horizontal_x1 > maxx) || (horizontal_x2 <= minx))
     168           0 :                         continue;
     169             : 
     170             :                     // Fill the horizontal segment (separately from the rest).
     171          36 :                     if (bAvoidBurningSamePoints)
     172             :                     {
     173          10 :                         polyInts2[ints2++] = horizontal_x1;
     174          10 :                         polyInts2[ints2++] = horizontal_x2;
     175             :                     }
     176             :                     else
     177             :                     {
     178          26 :                         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         116 :                 continue;
     186             :             }
     187             : 
     188       18846 :             if (dy < dy2 && dy >= dy1)
     189             :             {
     190       18766 :                 const double intersect =
     191       18766 :                     (dy - dy1) * (dx2 - dx1) / (dy2 - dy1) + dx1;
     192             : 
     193       18766 :                 polyInts[ints++] = static_cast<int>(floor(intersect + 0.5));
     194             :             }
     195             :         }
     196             : 
     197        9425 :         std::sort(polyInts.begin(), polyInts.begin() + ints);
     198        9425 :         std::sort(polyInts2.begin(), polyInts2.begin() + ints2);
     199             : 
     200       18808 :         for (int i = 0; i + 1 < ints; i += 2)
     201             :         {
     202        9383 :             if (polyInts[i] <= maxx && polyInts[i + 1] > minx)
     203             :             {
     204        8713 :                 pfnScanlineFunc(pCBData, y, polyInts[i], polyInts[i + 1] - 1,
     205             :                                 dfVariant == nullptr ? 0 : dfVariant[0]);
     206             :             }
     207             :         }
     208             : 
     209        9435 :         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          89 : 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          89 :     constexpr double EPSILON_INTERSECT_ONLY = 1e-4;
     395             : 
     396          89 :     if (!nPartCount)
     397           0 :         return;
     398             : 
     399         181 :     for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
     400             :     {
     401         184 :         std::set<std::pair<int, int>> lastBurntPoints;
     402         184 :         std::set<std::pair<int, int>> newBurntPoints;
     403             : 
     404         710 :         for (int j = 1; j < panPartSize[i]; j++)
     405             :         {
     406         618 :             lastBurntPoints = std::move(newBurntPoints);
     407         618 :             newBurntPoints.clear();
     408             : 
     409         618 :             double dfX = padfX[n + j - 1];
     410         618 :             double dfY = padfY[n + j - 1];
     411             : 
     412         618 :             double dfXEnd = padfX[n + j];
     413         618 :             double dfYEnd = padfY[n + j];
     414             : 
     415         618 :             double dfVariant = 0.0;
     416         618 :             double dfVariantEnd = 0.0;
     417         618 :             if (padfVariant != nullptr &&
     418          11 :                 static_cast<GDALRasterizeInfo *>(pCBData)->eBurnValueSource !=
     419             :                     GBV_UserBurnValue)
     420             :             {
     421          11 :                 dfVariant = padfVariant[n + j - 1];
     422          11 :                 dfVariantEnd = padfVariant[n + j];
     423             :             }
     424             : 
     425             :             // Skip segments that are off the target region.
     426         618 :             if ((dfY < 0.0 && dfYEnd < 0.0) ||
     427         616 :                 (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
     428         614 :                 (dfX < 0.0 && dfXEnd < 0.0) ||
     429         592 :                 (dfX > nRasterXSize && dfXEnd > nRasterXSize))
     430         461 :                 continue;
     431             : 
     432             :             // Swap if needed so we can proceed from left2right (X increasing)
     433         573 :             if (dfX > dfXEnd)
     434             :             {
     435         212 :                 std::swap(dfX, dfXEnd);
     436         212 :                 std::swap(dfY, dfYEnd);
     437         212 :                 std::swap(dfVariant, dfVariantEnd);
     438             :             }
     439             : 
     440             :             // Special case for vertical lines.
     441             : 
     442         573 :             if (fabs(dfX - dfXEnd) < .01)
     443             :             {
     444         189 :                 if (bIntersectOnly)
     445             :                 {
     446         180 :                     if (std::abs(dfX - std::round(dfX)) <
     447         252 :                             EPSILON_INTERSECT_ONLY &&
     448          72 :                         std::abs(dfXEnd - std::round(dfXEnd)) <
     449             :                             EPSILON_INTERSECT_ONLY)
     450          71 :                         continue;
     451             :                 }
     452             : 
     453         118 :                 if (dfYEnd < dfY)
     454             :                 {
     455          61 :                     std::swap(dfY, dfYEnd);
     456          61 :                     std::swap(dfVariant, dfVariantEnd);
     457             :                 }
     458             : 
     459         118 :                 const int iX = static_cast<int>(floor(dfXEnd));
     460         118 :                 int iY = static_cast<int>(floor(dfY));
     461         118 :                 int iYEnd =
     462         118 :                     static_cast<int>(floor(dfYEnd - EPSILON_INTERSECT_ONLY));
     463             : 
     464         118 :                 if (iX < 0 || iX >= nRasterXSize)
     465           1 :                     continue;
     466             : 
     467         117 :                 double dfDeltaVariant = 0.0;
     468         117 :                 if (dfYEnd - dfY > 0.0)
     469         117 :                     dfDeltaVariant = (dfVariantEnd - dfVariant) /
     470         117 :                                      (dfYEnd - dfY);  // Per unit change in iY.
     471             : 
     472             :                 // Clip to the borders of the target region.
     473         117 :                 if (iY < 0)
     474           0 :                     iY = 0;
     475         117 :                 if (iYEnd >= nRasterYSize)
     476           0 :                     iYEnd = nRasterYSize - 1;
     477         117 :                 dfVariant += dfDeltaVariant * (iY - dfY);
     478             : 
     479         117 :                 if (padfVariant == nullptr)
     480             :                 {
     481        1203 :                     for (; iY <= iYEnd; iY++)
     482             :                     {
     483        1090 :                         if (bAvoidBurningSamePoints)
     484             :                         {
     485         146 :                             auto yx = std::pair<int, int>(iY, iX);
     486         146 :                             if (lastBurntPoints.find(yx) !=
     487         292 :                                 lastBurntPoints.end())
     488             :                             {
     489          13 :                                 continue;
     490             :                             }
     491         133 :                             newBurntPoints.insert(yx);
     492             :                         }
     493        1077 :                         pfnPointFunc(pCBData, iY, iX, 0.0);
     494             :                     }
     495             :                 }
     496             :                 else
     497             :                 {
     498          14 :                     for (; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant)
     499             :                     {
     500          10 :                         if (bAvoidBurningSamePoints)
     501             :                         {
     502           0 :                             auto yx = std::pair<int, int>(iY, iX);
     503           0 :                             if (lastBurntPoints.find(yx) !=
     504           0 :                                 lastBurntPoints.end())
     505             :                             {
     506           0 :                                 continue;
     507             :                             }
     508           0 :                             newBurntPoints.insert(yx);
     509             :                         }
     510          10 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     511             :                     }
     512             :                 }
     513             : 
     514         117 :                 continue;  // Next segment.
     515             :             }
     516             : 
     517         384 :             const double dfDeltaVariant =
     518         384 :                 (dfVariantEnd - dfVariant) /
     519         384 :                 (dfXEnd - dfX);  // Per unit change in iX.
     520             : 
     521             :             // Special case for horizontal lines.
     522         384 :             if (fabs(dfY - dfYEnd) < .01)
     523             :             {
     524         227 :                 if (bIntersectOnly)
     525             :                 {
     526         218 :                     if (std::abs(dfY - std::round(dfY)) <
     527         288 :                             EPSILON_INTERSECT_ONLY &&
     528          70 :                         std::abs(dfYEnd - std::round(dfYEnd)) <
     529             :                             EPSILON_INTERSECT_ONLY)
     530          69 :                         continue;
     531             :                 }
     532             : 
     533         158 :                 if (dfXEnd < dfX)
     534             :                 {
     535           0 :                     std::swap(dfX, dfXEnd);
     536           0 :                     std::swap(dfVariant, dfVariantEnd);
     537             :                 }
     538             : 
     539         158 :                 int iX = static_cast<int>(floor(dfX));
     540         158 :                 const int iY = static_cast<int>(floor(dfY));
     541         158 :                 int iXEnd =
     542         158 :                     static_cast<int>(floor(dfXEnd - EPSILON_INTERSECT_ONLY));
     543             : 
     544         158 :                 if (iY < 0 || iY >= nRasterYSize)
     545           0 :                     continue;
     546             : 
     547             :                 // Clip to the borders of the target region.
     548         158 :                 if (iX < 0)
     549           2 :                     iX = 0;
     550         158 :                 if (iXEnd >= nRasterXSize)
     551           1 :                     iXEnd = nRasterXSize - 1;
     552         158 :                 dfVariant += dfDeltaVariant * (iX - dfX);
     553             : 
     554         158 :                 if (padfVariant == nullptr)
     555             :                 {
     556        1086 :                     for (; iX <= iXEnd; iX++)
     557             :                     {
     558         932 :                         if (bAvoidBurningSamePoints)
     559             :                         {
     560         144 :                             auto yx = std::pair<int, int>(iY, iX);
     561         144 :                             if (lastBurntPoints.find(yx) !=
     562         288 :                                 lastBurntPoints.end())
     563             :                             {
     564          12 :                                 continue;
     565             :                             }
     566         132 :                             newBurntPoints.insert(yx);
     567             :                         }
     568         920 :                         pfnPointFunc(pCBData, iY, iX, 0.0);
     569             :                     }
     570             :                 }
     571             :                 else
     572             :                 {
     573          14 :                     for (; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant)
     574             :                     {
     575          10 :                         if (bAvoidBurningSamePoints)
     576             :                         {
     577           0 :                             auto yx = std::pair<int, int>(iY, iX);
     578           0 :                             if (lastBurntPoints.find(yx) !=
     579           0 :                                 lastBurntPoints.end())
     580             :                             {
     581           0 :                                 continue;
     582             :                             }
     583           0 :                             newBurntPoints.insert(yx);
     584             :                         }
     585          10 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     586             :                     }
     587             :                 }
     588             : 
     589         158 :                 continue;  // Next segment.
     590             :             }
     591             : 
     592             :             /* --------------------------------------------------------------------
     593             :              */
     594             :             /*      General case - left to right sloped. */
     595             :             /* --------------------------------------------------------------------
     596             :              */
     597         157 :             const double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
     598             : 
     599             :             // Clip segment in X.
     600         157 :             if (dfXEnd > nRasterXSize)
     601             :             {
     602           0 :                 dfYEnd -= (dfXEnd - nRasterXSize) * dfSlope;
     603           0 :                 dfXEnd = nRasterXSize;
     604             :             }
     605         157 :             if (dfX < 0.0)
     606             :             {
     607           2 :                 dfY += (0.0 - dfX) * dfSlope;
     608           2 :                 dfVariant += dfDeltaVariant * (0.0 - dfX);
     609           2 :                 dfX = 0.0;
     610             :             }
     611             : 
     612             :             // Clip segment in Y.
     613         157 :             if (dfYEnd > dfY)
     614             :             {
     615           9 :                 if (dfY < 0.0)
     616             :                 {
     617           0 :                     const double dfDiffX = (0.0 - dfY) / dfSlope;
     618           0 :                     dfX += dfDiffX;
     619           0 :                     dfVariant += dfDeltaVariant * dfDiffX;
     620           0 :                     dfY = 0.0;
     621             :                 }
     622           9 :                 if (dfYEnd >= nRasterYSize)
     623             :                 {
     624           2 :                     dfXEnd += (dfYEnd - nRasterYSize) / dfSlope;
     625           2 :                     if (dfXEnd > nRasterXSize)
     626           2 :                         dfXEnd = nRasterXSize;
     627             :                     // dfYEnd is no longer used afterwards, but for
     628             :                     // consistency it should be:
     629             :                     // dfYEnd = nRasterXSize;
     630             :                 }
     631             :             }
     632             :             else
     633             :             {
     634         148 :                 if (dfY >= nRasterYSize)
     635             :                 {
     636           2 :                     const double dfDiffX = (nRasterYSize - dfY) / dfSlope;
     637           2 :                     dfX += dfDiffX;
     638           2 :                     dfVariant += dfDeltaVariant * dfDiffX;
     639           2 :                     dfY = nRasterYSize;
     640             :                 }
     641         148 :                 if (dfYEnd < 0.0)
     642             :                 {
     643           0 :                     dfXEnd -= (dfYEnd - 0) / dfSlope;
     644             :                     // dfYEnd is no longer used afterwards, but for
     645             :                     // consistency it should be:
     646             :                     // dfYEnd = 0.0;
     647             :                 }
     648             :             }
     649             : 
     650             :             // Step from pixel to pixel.
     651        5157 :             while (dfX >= 0.0 && dfX < dfXEnd)
     652             :             {
     653        5000 :                 const int iX = static_cast<int>(floor(dfX));
     654        5000 :                 const int iY = static_cast<int>(floor(dfY));
     655             : 
     656             :                 // Burn in the current point.
     657             :                 // We should be able to drop the Y check because we clipped
     658             :                 // in Y, but there may be some error with all the small steps.
     659        5000 :                 if (iY >= 0 && iY < nRasterYSize)
     660             :                 {
     661        4997 :                     if (bAvoidBurningSamePoints)
     662             :                     {
     663          42 :                         auto yx = std::pair<int, int>(iY, iX);
     664          79 :                         if (lastBurntPoints.find(yx) == lastBurntPoints.end() &&
     665          79 :                             newBurntPoints.find(yx) == newBurntPoints.end())
     666             :                         {
     667          27 :                             newBurntPoints.insert(yx);
     668          27 :                             pfnPointFunc(pCBData, iY, iX, dfVariant);
     669             :                         }
     670             :                     }
     671             :                     else
     672             :                     {
     673        4955 :                         pfnPointFunc(pCBData, iY, iX, dfVariant);
     674             :                     }
     675             :                 }
     676             : 
     677        5000 :                 double dfStepX = floor(dfX + 1.0) - dfX;
     678        5000 :                 double dfStepY = dfStepX * dfSlope;
     679             : 
     680             :                 // Step to right pixel without changing scanline?
     681        5000 :                 if (static_cast<int>(floor(dfY + dfStepY)) == iY)
     682             :                 {
     683        3149 :                     dfX += dfStepX;
     684        3149 :                     dfY += dfStepY;
     685        3149 :                     dfVariant += dfDeltaVariant * dfStepX;
     686             :                 }
     687        1851 :                 else if (dfSlope < 0)
     688             :                 {
     689        1306 :                     dfStepY = iY - dfY;
     690        1306 :                     if (dfStepY > -0.000000001)
     691         651 :                         dfStepY = -0.000000001;
     692             : 
     693        1306 :                     dfStepX = dfStepY / dfSlope;
     694        1306 :                     dfX += dfStepX;
     695        1306 :                     dfY += dfStepY;
     696        1306 :                     dfVariant += dfDeltaVariant * dfStepX;
     697             :                 }
     698             :                 else
     699             :                 {
     700         545 :                     dfStepY = (iY + 1) - dfY;
     701         545 :                     if (dfStepY < 0.000000001)
     702           1 :                         dfStepY = 0.000000001;
     703             : 
     704         545 :                     dfStepX = dfStepY / dfSlope;
     705         545 :                     dfX += dfStepX;
     706         545 :                     dfY += dfStepY;
     707         545 :                     dfVariant += dfDeltaVariant * dfStepX;
     708             :                 }
     709             :             }  // Next step along segment.
     710             :         }      // Next segment.
     711             :     }          // Next part.
     712             : }

Generated by: LCOV version 1.14