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

Generated by: LCOV version 1.14