/src/gdal/alg/llrasterize.cpp

Source
/******************************************************************************
 *
 * Project:  GDAL
 * Purpose:  Vector polygon rasterization code.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 2000, Frank Warmerdam <warmerdam@pobox.com>
 * Copyright (c) 2011, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_port.h"
#include "gdal_alg_priv.h"

#include <cmath>
#include <cstdlib>
#include <cstring>

#include <algorithm>
#include <set>
#include <utility>
#include <vector>

#include "gdal_alg.h"

/************************************************************************/
/*                       dllImageFilledPolygon()                        */
/*                                                                      */
/*      Perform scanline conversion of the passed multi-ring            */
/*      polygon.  Note the polygon does not need to be explicitly       */
/*      closed.  The scanline function will be called with              */
/*      horizontal scanline chunks which may not be entirely            */
/*      contained within the valid raster area (in the X                */
/*      direction).                                                     */
/*                                                                      */
/*      NEW: Nodes' coordinate are kept as double  in order             */
/*      to compute accurately the intersections with the lines          */
/*                                                                      */
/*      A pixel is considered inside a polygon if its center            */
/*      falls inside the polygon. This is robust unless                 */
/*      the nodes are placed in the center of the pixels in which       */
/*      case, due to numerical inaccuracies, it's hard to predict       */
/*      if the pixel will be considered inside or outside the shape.    */
/************************************************************************/

/*
 * NOTE: This code was originally adapted from the gdImageFilledPolygon()
 * function in libgd.
 *
 * http://www.boutell.com/gd/
 *
 * It was later adapted for direct inclusion in GDAL and relicensed under
 * the GDAL MIT license (pulled from the OpenEV distribution).
 */

void GDALdllImageFilledPolygon(int nRasterXSize, int nRasterYSize,
                               int nPartCount, const int *panPartSize,
                               const double *padfX, const double *padfY,
                               const double *dfVariant,
                               llScanlineFunc pfnScanlineFunc,
                               GDALRasterizeInfo *pCBData,
                               bool bAvoidBurningSamePoints)
{
    if (!nPartCount)
    {
        return;
    }

    int n = 0;
    for (int part = 0; part < nPartCount; part++)
        n += panPartSize[part];

    std::vector<int> polyInts(n);
    std::vector<int> polyInts2;
    if (bAvoidBurningSamePoints)
        polyInts2.resize(n);

    double dminy = padfY[0];
    double dmaxy = padfY[0];
    for (int i = 1; i < n; i++)
    {
        if (padfY[i] < dminy)
        {
            dminy = padfY[i];
        }
        else if (padfY[i] > dmaxy)
        {
            dmaxy = padfY[i];
        }
    }
    const int miny = static_cast<int>(std::max(0.0, dminy));
    const int maxy =
        static_cast<int>(std::min<double>(dmaxy, nRasterYSize - 1));

    constexpr int minx = 0;
    const int maxx = nRasterXSize - 1;

    // Fix in 1.3: count a vertex only once.
    for (int y = miny; y <= maxy; y++)
    {
        int partoffset = 0;

        const double dy = y + 0.5;  // Center height of line.

        int part = 0;
        int ints = 0;
        int ints2 = 0;

        for (int i = 0; i < n; i++)
        {
            if (i == partoffset + panPartSize[part])
            {
                partoffset += panPartSize[part];
                part++;
            }

            int ind1 = 0;
            int ind2 = 0;
            if (i == partoffset)
            {
                ind1 = partoffset + panPartSize[part] - 1;
                ind2 = partoffset;
            }
            else
            {
                ind1 = i - 1;
                ind2 = i;
            }

            double dy1 = padfY[ind1];
            double dy2 = padfY[ind2];

            if ((dy1 < dy && dy2 < dy) || (dy1 > dy && dy2 > dy))
                continue;

            double dx1 = 0.0;
            double dx2 = 0.0;
            if (dy1 < dy2)
            {
                dx1 = padfX[ind1];
                dx2 = padfX[ind2];
            }
            else if (dy1 > dy2)
            {
                std::swap(dy1, dy2);
                dx2 = padfX[ind1];
                dx1 = padfX[ind2];
            }
            else  // if( fabs(dy1-dy2) < 1.e-6 )
            {

                // AE: DO NOT skip bottom horizontal segments
                // -Fill them separately-
                if (padfX[ind1] > padfX[ind2])
                {
                    const double dhorizontal_x1 = floor(padfX[ind2] + 0.5);
                    const double dhorizontal_x2 = floor(padfX[ind1] + 0.5);

                    if ((dhorizontal_x1 > static_cast<double>(maxx)) ||
                        (dhorizontal_x2 <= 0))
                        continue;
                    const int horizontal_x1 =
                        static_cast<int>(std::max(dhorizontal_x1, 0.0));
                    const int horizontal_x2 = static_cast<int>(
                        std::min<double>(dhorizontal_x2, nRasterXSize));

                    // Fill the horizontal segment (separately from the rest).
                    if (bAvoidBurningSamePoints)
                    {
                        polyInts2[ints2++] = horizontal_x1;
                        polyInts2[ints2++] = horizontal_x2;
                    }
                    else
                    {
                        pfnScanlineFunc(
                            pCBData, y, horizontal_x1, horizontal_x2 - 1,
                            dfVariant == nullptr ? 0 : dfVariant[0]);
                    }
                }
                // else: Skip top horizontal segments.
                // They are already filled in the regular loop.
                continue;
            }

            if (dy < dy2 && dy >= dy1)
            {
                const double intersect = std::clamp<double>(
                    (dy - dy1) * (dx2 - dx1) / (dy2 - dy1) + dx1, INT_MIN,
                    INT_MAX);

                polyInts[ints++] = static_cast<int>(floor(intersect + 0.5));
            }
        }

        std::sort(polyInts.begin(), polyInts.begin() + ints);
        std::sort(polyInts2.begin(), polyInts2.begin() + ints2);

        for (int i = 0; i + 1 < ints; i += 2)
        {
            if (polyInts[i] <= maxx && polyInts[i + 1] > minx)
            {
                pfnScanlineFunc(pCBData, y, polyInts[i], polyInts[i + 1] - 1,
                                dfVariant == nullptr ? 0 : dfVariant[0]);
            }
        }

        for (int i2 = 0, i = 0; i2 + 1 < ints2; i2 += 2)
        {
            if (polyInts2[i2] <= maxx && polyInts2[i2 + 1] > minx)
            {
                // "synchronize" polyInts[i] with polyInts2[i2]
                while (i + 1 < ints && polyInts[i] < polyInts2[i2])
                    i += 2;
                // Only burn if we don't have a common segment between
                // polyInts[] and polyInts2[]
                if (i + 1 >= ints || polyInts[i] != polyInts2[i2])
                {
                    pfnScanlineFunc(pCBData, y, polyInts2[i2],
                                    polyInts2[i2 + 1] - 1,
                                    dfVariant == nullptr ? 0 : dfVariant[0]);
                }
            }
        }
    }
}

/************************************************************************/
/*                         GDALdllImagePoint()                          */
/************************************************************************/

void GDALdllImagePoint(int nRasterXSize, int nRasterYSize, int nPartCount,
                       const int * /*panPartSize*/, const double *padfX,
                       const double *padfY, const double *padfVariant,
                       llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData)
{
    for (int i = 0; i < nPartCount; i++)
    {
        const double dfX = padfX[i];
        const double dfY = padfY[i];
        double dfVariant = 0.0;
        if (padfVariant != nullptr)
            dfVariant = padfVariant[i];

        if (0 <= dfX && dfX < nRasterXSize && 0 <= dfY && dfY < nRasterYSize)
            pfnPointFunc(pCBData, static_cast<int>(dfY), static_cast<int>(dfX),
                         dfVariant);
    }
}

/************************************************************************/
/*                         GDALdllImageLine()                           */
/************************************************************************/

void GDALdllImageLine(int nRasterXSize, int nRasterYSize, int nPartCount,
                      const int *panPartSize, const double *padfX,
                      const double *padfY, const double *padfVariant,
                      llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData)
{
    if (!nPartCount)
        return;

    for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
    {
        for (int j = 1; j < panPartSize[i]; j++)
        {
            double dfX = padfX[n + j - 1];
            double dfY = padfY[n + j - 1];
            double dfXEnd = padfX[n + j];
            double dfYEnd = padfY[n + j];

            // Skip segments that are off the target region.
            if ((dfY < 0.0 && dfYEnd < 0.0) ||
                (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
                (dfX < 0.0 && dfXEnd < 0.0) ||
                (dfX > nRasterXSize && dfXEnd > nRasterXSize))
                continue;

            // TODO: clamp coordinates to [0, nRasterXSize] * [0, nRasterYSize]
            if (!(dfX >= INT_MIN && dfX <= INT_MAX && dfY >= INT_MIN &&
                  dfY <= INT_MAX && dfXEnd >= INT_MIN && dfXEnd <= INT_MAX &&
                  dfYEnd >= INT_MIN && dfYEnd <= INT_MAX))
            {
                CPLErrorOnce(
                    CE_Warning, CPLE_AppDefined,
                    "GDALdllImageLine(): coordinates outside of int range");
                continue;
            }

            int iX = static_cast<int>(floor(dfX));
            int iY = static_cast<int>(floor(dfY));

            const int iX1 = static_cast<int>(floor(dfXEnd));
            const int iY1 = static_cast<int>(floor(dfYEnd));

            double dfVariant = 0.0;
            double dfVariant1 = 0.0;
            if (padfVariant != nullptr &&
                pCBData->eBurnValueSource != GBV_UserBurnValue)
            {
                dfVariant = padfVariant[n + j - 1];
                dfVariant1 = padfVariant[n + j];
            }

            int nDeltaX = std::abs(iX1 - iX);
            int nDeltaY = std::abs(iY1 - iY);

            // Step direction depends on line direction.
            const int nXStep = (iX > iX1) ? -1 : 1;
            const int nYStep = (iY > iY1) ? -1 : 1;

            // Determine the line slope.
            if (nDeltaX >= nDeltaY)
            {
                const int nXError = nDeltaY << 1;
                const int nYError = nXError - (nDeltaX << 1);
                int nError = nXError - nDeltaX;
                // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
                // but if it is == 0, dfDeltaVariant is not really used, so any
                // value is okay.
                const double dfDeltaVariant =
                    nDeltaX == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaX;

                // Do not burn the end point, unless we are in the last
                // segment. This is to avoid burning twice intermediate points,
                // which causes artifacts in Add mode
                if (j != panPartSize[i] - 1)
                {
                    nDeltaX--;
                }

                while (nDeltaX-- >= 0)
                {
                    if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
                        iY < nRasterYSize)
                        pfnPointFunc(pCBData, iY, iX, dfVariant);

                    dfVariant += dfDeltaVariant;
                    iX += nXStep;
                    if (nError > 0)
                    {
                        iY += nYStep;
                        nError += nYError;
                    }
                    else
                    {
                        nError += nXError;
                    }
                }
            }
            else
            {
                const int nXError = nDeltaX << 1;
                const int nYError = nXError - (nDeltaY << 1);
                int nError = nXError - nDeltaY;
                // == 0 makes clang -fcatch-undefined-behavior -ftrapv happy,
                // but if it is == 0, dfDeltaVariant is not really used, so any
                // value is okay.
                double dfDeltaVariant =
                    nDeltaY == 0 ? 0.0 : (dfVariant1 - dfVariant) / nDeltaY;

                // Do not burn the end point, unless we are in the last
                // segment. This is to avoid burning twice intermediate points,
                // which causes artifacts in Add mode
                if (j != panPartSize[i] - 1)
                {
                    nDeltaY--;
                }

                while (nDeltaY-- >= 0)
                {
                    if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
                        iY < nRasterYSize)
                        pfnPointFunc(pCBData, iY, iX, dfVariant);

                    dfVariant += dfDeltaVariant;
                    iY += nYStep;
                    if (nError > 0)
                    {
                        iX += nXStep;
                        nError += nYError;
                    }
                    else
                    {
                        nError += nXError;
                    }
                }
            }
        }
    }
}

/************************************************************************/
/*                     GDALdllImageLineAllTouched()                     */
/*                                                                      */
/*      This alternate line drawing algorithm attempts to ensure        */
/*      that every pixel touched at all by the line will get set.       */
/*      @param padfVariant should contain the values that are to be     */
/*      added to the burn value.  The values along the line between the */
/*      points will be linearly interpolated. These values are used     */
/*      only if pCBData->eBurnValueSource is set to something other     */
/*      than GBV_UserBurnValue. If NULL is passed, a monotonous line    */
/*      will be drawn with the burn value.                              */
/************************************************************************/

void GDALdllImageLineAllTouched(
    int nRasterXSize, int nRasterYSize, int nPartCount, const int *panPartSize,
    const double *padfX, const double *padfY, const double *padfVariant,
    llPointFunc pfnPointFunc, GDALRasterizeInfo *pCBData,
    bool bAvoidBurningSamePoints, bool bIntersectOnly)

{
    // This is an epsilon to detect geometries that are aligned with pixel
    // coordinates. Hard to find the right value. We put it to that value
    // to satisfy the scenarios of https://github.com/OSGeo/gdal/issues/7523
    // and https://github.com/OSGeo/gdal/issues/6414
    constexpr double EPSILON_INTERSECT_ONLY = 1e-4;

    if (!nPartCount)
        return;

    for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
    {
        std::set<std::pair<int, int>> lastBurntPoints;
        std::set<std::pair<int, int>> newBurntPoints;

        for (int j = 1; j < panPartSize[i]; j++)
        {
            lastBurntPoints = std::move(newBurntPoints);
            newBurntPoints.clear();

            double dfX = padfX[n + j - 1];
            double dfY = padfY[n + j - 1];

            double dfXEnd = padfX[n + j];
            double dfYEnd = padfY[n + j];

            // Skip segments that are off the target region.
            if ((dfY < 0.0 && dfYEnd < 0.0) ||
                (dfY > nRasterYSize && dfYEnd > nRasterYSize) ||
                (dfX < 0.0 && dfXEnd < 0.0) ||
                (dfX > nRasterXSize && dfXEnd > nRasterXSize))
                continue;

            // TODO: clamp coordinates to [0, nRasterXSize] * [0, nRasterYSize]
            if (!(dfX >= INT_MIN && dfX <= INT_MAX && dfY >= INT_MIN &&
                  dfY <= INT_MAX && dfXEnd >= INT_MIN && dfXEnd <= INT_MAX &&
                  dfYEnd >= INT_MIN && dfYEnd <= INT_MAX))
            {
                CPLErrorOnce(CE_Warning, CPLE_AppDefined,
                             "GDALdllImageLineAllTouched(): coordinates "
                             "outside of int range");
                continue;
            }

            double dfVariant = 0.0;
            double dfVariantEnd = 0.0;
            if (padfVariant != nullptr &&
                static_cast<GDALRasterizeInfo *>(pCBData)->eBurnValueSource !=
                    GBV_UserBurnValue)
            {
                dfVariant = padfVariant[n + j - 1];
                dfVariantEnd = padfVariant[n + j];
            }

            // Swap if needed so we can proceed from left2right (X increasing)
            if (dfX > dfXEnd)
            {
                std::swap(dfX, dfXEnd);
                std::swap(dfY, dfYEnd);
                std::swap(dfVariant, dfVariantEnd);
            }

            // Special case for vertical lines.

            if (fabs(dfX - dfXEnd) < .01)
            {
                if (bIntersectOnly)
                {
                    if (std::abs(dfX - std::round(dfX)) <
                            EPSILON_INTERSECT_ONLY &&
                        std::abs(dfXEnd - std::round(dfXEnd)) <
                            EPSILON_INTERSECT_ONLY)
                        continue;
                }

                if (dfYEnd < dfY)
                {
                    std::swap(dfY, dfYEnd);
                    std::swap(dfVariant, dfVariantEnd);
                }

                const int iX = static_cast<int>(floor(dfXEnd));
                int iY = static_cast<int>(floor(dfY));
                int iYEnd =
                    static_cast<int>(floor(dfYEnd - EPSILON_INTERSECT_ONLY));

                if (iX < 0 || iX >= nRasterXSize)
                    continue;

                double dfDeltaVariant = 0.0;
                if (dfYEnd - dfY > 0.0)
                    dfDeltaVariant = (dfVariantEnd - dfVariant) /
                                     (dfYEnd - dfY);  // Per unit change in iY.

                // Clip to the borders of the target region.
                if (iY < 0)
                    iY = 0;
                if (iYEnd >= nRasterYSize)
                    iYEnd = nRasterYSize - 1;
                dfVariant += dfDeltaVariant * (iY - dfY);

                if (padfVariant == nullptr)
                {
                    for (; iY <= iYEnd; iY++)
                    {
                        if (bAvoidBurningSamePoints)
                        {
                            auto yx = std::pair<int, int>(iY, iX);
                            if (lastBurntPoints.find(yx) !=
                                lastBurntPoints.end())
                            {
                                continue;
                            }
                            newBurntPoints.insert(yx);
                        }
                        pfnPointFunc(pCBData, iY, iX, 0.0);
                    }
                }
                else
                {
                    for (; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant)
                    {
                        if (bAvoidBurningSamePoints)
                        {
                            auto yx = std::pair<int, int>(iY, iX);
                            if (lastBurntPoints.find(yx) !=
                                lastBurntPoints.end())
                            {
                                continue;
                            }
                            newBurntPoints.insert(yx);
                        }
                        pfnPointFunc(pCBData, iY, iX, dfVariant);
                    }
                }

                continue;  // Next segment.
            }

            const double dfDeltaVariant =
                (dfVariantEnd - dfVariant) /
                (dfXEnd - dfX);  // Per unit change in iX.

            // Special case for horizontal lines.
            if (fabs(dfY - dfYEnd) < .01)
            {
                if (bIntersectOnly)
                {
                    if (std::abs(dfY - std::round(dfY)) <
                            EPSILON_INTERSECT_ONLY &&
                        std::abs(dfYEnd - std::round(dfYEnd)) <
                            EPSILON_INTERSECT_ONLY)
                        continue;
                }

                if (dfXEnd < dfX)
                {
                    std::swap(dfX, dfXEnd);
                    std::swap(dfVariant, dfVariantEnd);
                }

                int iX = static_cast<int>(floor(dfX));
                const int iY = static_cast<int>(floor(dfY));
                int iXEnd =
                    static_cast<int>(floor(dfXEnd - EPSILON_INTERSECT_ONLY));

                if (iY < 0 || iY >= nRasterYSize)
                    continue;

                // Clip to the borders of the target region.
                if (iX < 0)
                    iX = 0;
                if (iXEnd >= nRasterXSize)
                    iXEnd = nRasterXSize - 1;
                dfVariant += dfDeltaVariant * (iX - dfX);

                if (padfVariant == nullptr)
                {
                    for (; iX <= iXEnd; iX++)
                    {
                        if (bAvoidBurningSamePoints)
                        {
                            auto yx = std::pair<int, int>(iY, iX);
                            if (lastBurntPoints.find(yx) !=
                                lastBurntPoints.end())
                            {
                                continue;
                            }
                            newBurntPoints.insert(yx);
                        }
                        pfnPointFunc(pCBData, iY, iX, 0.0);
                    }
                }
                else
                {
                    for (; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant)
                    {
                        if (bAvoidBurningSamePoints)
                        {
                            auto yx = std::pair<int, int>(iY, iX);
                            if (lastBurntPoints.find(yx) !=
                                lastBurntPoints.end())
                            {
                                continue;
                            }
                            newBurntPoints.insert(yx);
                        }
                        pfnPointFunc(pCBData, iY, iX, dfVariant);
                    }
                }

                continue;  // Next segment.
            }

            /* --------------------------------------------------------------------
             */
            /*      General case - left to right sloped. */
            /* --------------------------------------------------------------------
             */
            const double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);

            // Clip segment in X.
            if (dfXEnd > nRasterXSize)
            {
                dfYEnd -= (dfXEnd - nRasterXSize) * dfSlope;
                dfXEnd = nRasterXSize;
            }
            if (dfX < 0.0)
            {
                dfY += (0.0 - dfX) * dfSlope;
                dfVariant += dfDeltaVariant * (0.0 - dfX);
                dfX = 0.0;
            }

            // Clip segment in Y.
            if (dfYEnd > dfY)
            {
                if (dfY < 0.0)
                {
                    const double dfDiffX = (0.0 - dfY) / dfSlope;
                    dfX += dfDiffX;
                    dfVariant += dfDeltaVariant * dfDiffX;
                    dfY = 0.0;
                }
                if (dfYEnd >= nRasterYSize)
                {
                    dfXEnd += (dfYEnd - nRasterYSize) / dfSlope;
                    if (dfXEnd > nRasterXSize)
                        dfXEnd = nRasterXSize;
                    // dfYEnd is no longer used afterwards, but for
                    // consistency it should be:
                    // dfYEnd = nRasterXSize;
                }
            }
            else
            {
                if (dfY >= nRasterYSize)
                {
                    const double dfDiffX = (nRasterYSize - dfY) / dfSlope;
                    dfX += dfDiffX;
                    dfVariant += dfDeltaVariant * dfDiffX;
                    dfY = nRasterYSize;
                }
                if (dfYEnd < 0.0)
                {
                    dfXEnd -= (dfYEnd - 0) / dfSlope;
                    // dfYEnd is no longer used afterwards, but for
                    // consistency it should be:
                    // dfYEnd = 0.0;
                }
            }

            // Step from pixel to pixel.
            while (dfX >= 0.0 && dfX < dfXEnd)
            {
                const int iX = static_cast<int>(floor(dfX));
                const int iY = static_cast<int>(floor(dfY));

                // Burn in the current point.
                // We should be able to drop the Y check because we clipped
                // in Y, but there may be some error with all the small steps.
                if (iY >= 0 && iY < nRasterYSize)
                {
                    if (bAvoidBurningSamePoints)
                    {
                        auto yx = std::pair<int, int>(iY, iX);
                        if (lastBurntPoints.find(yx) == lastBurntPoints.end() &&
                            newBurntPoints.find(yx) == newBurntPoints.end())
                        {
                            newBurntPoints.insert(yx);
                            pfnPointFunc(pCBData, iY, iX, dfVariant);
                        }
                    }
                    else
                    {
                        pfnPointFunc(pCBData, iY, iX, dfVariant);
                    }
                }

                double dfStepX = floor(dfX + 1.0) - dfX;
                double dfStepY = dfStepX * dfSlope;

                // Step to right pixel without changing scanline?
                if (static_cast<int>(floor(dfY + dfStepY)) == iY)
                {
                    dfX += dfStepX;
                    dfY += dfStepY;
                    dfVariant += dfDeltaVariant * dfStepX;
                }
                else if (dfSlope < 0)
                {
                    dfStepY = iY - dfY;
                    if (dfStepY > -0.000000001)
                        dfStepY = -0.000000001;

                    dfStepX = dfStepY / dfSlope;
                    dfX += dfStepX;
                    dfY += dfStepY;
                    dfVariant += dfDeltaVariant * dfStepX;
                }
                else
                {
                    dfStepY = (iY + 1) - dfY;
                    if (dfStepY < 0.000000001)
                        dfStepY = 0.000000001;

                    dfStepX = dfStepY / dfSlope;
                    dfX += dfStepX;
                    dfY += dfStepY;
                    dfVariant += dfDeltaVariant * dfStepX;
                }
            }  // Next step along segment.
        }      // Next segment.
    }          // Next part.
}

Coverage Report

Created: 2025-11-16 06:25

Line	Count	Source
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		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	0	{
66	0	if (!nPartCount)
67	0	{
68	0	return;
69	0	}
70
71	0	int n = 0;
72	0	for (int part = 0; part < nPartCount; part++)
73	0	n += panPartSize[part];
74
75	0	std::vector<int> polyInts(n);
76	0	std::vector<int> polyInts2;
77	0	if (bAvoidBurningSamePoints)
78	0	polyInts2.resize(n);
79
80	0	double dminy = padfY[0];
81	0	double dmaxy = padfY[0];
82	0	for (int i = 1; i < n; i++)
83	0	{
84	0	if (padfY[i] < dminy)
85	0	{
86	0	dminy = padfY[i];
87	0	}
88	0	else if (padfY[i] > dmaxy)
89	0	{
90	0	dmaxy = padfY[i];
91	0	}
92	0	}
93	0	const int miny = static_cast<int>(std::max(0.0, dminy));
94	0	const int maxy =
95	0	static_cast<int>(std::min<double>(dmaxy, nRasterYSize - 1));
96
97	0	constexpr int minx = 0;
98	0	const int maxx = nRasterXSize - 1;
99
100		// Fix in 1.3: count a vertex only once.
101	0	for (int y = miny; y <= maxy; y++)
102	0	{
103	0	int partoffset = 0;
104
105	0	const double dy = y + 0.5; // Center height of line.
106
107	0	int part = 0;
108	0	int ints = 0;
109	0	int ints2 = 0;
110
111	0	for (int i = 0; i < n; i++)
112	0	{
113	0	if (i == partoffset + panPartSize[part])
114	0	{
115	0	partoffset += panPartSize[part];
116	0	part++;
117	0	}
118
119	0	int ind1 = 0;
120	0	int ind2 = 0;
121	0	if (i == partoffset)
122	0	{
123	0	ind1 = partoffset + panPartSize[part] - 1;
124	0	ind2 = partoffset;
125	0	}
126	0	else
127	0	{
128	0	ind1 = i - 1;
129	0	ind2 = i;
130	0	}
131
132	0	double dy1 = padfY[ind1];
133	0	double dy2 = padfY[ind2];
134
135	0	if ((dy1 < dy && dy2 < dy) \|\| (dy1 > dy && dy2 > dy))
136	0	continue;
137
138	0	double dx1 = 0.0;
139	0	double dx2 = 0.0;
140	0	if (dy1 < dy2)
141	0	{
142	0	dx1 = padfX[ind1];
143	0	dx2 = padfX[ind2];
144	0	}
145	0	else if (dy1 > dy2)
146	0	{
147	0	std::swap(dy1, dy2);
148	0	dx2 = padfX[ind1];
149	0	dx1 = padfX[ind2];
150	0	}
151	0	else // if( fabs(dy1-dy2) < 1.e-6 )
152	0	{
153
154		// AE: DO NOT skip bottom horizontal segments
155		// -Fill them separately-
156	0	if (padfX[ind1] > padfX[ind2])
157	0	{
158	0	const double dhorizontal_x1 = floor(padfX[ind2] + 0.5);
159	0	const double dhorizontal_x2 = floor(padfX[ind1] + 0.5);
160
161	0	if ((dhorizontal_x1 > static_cast<double>(maxx)) \|\|
162	0	(dhorizontal_x2 <= 0))
163	0	continue;
164	0	const int horizontal_x1 =
165	0	static_cast<int>(std::max(dhorizontal_x1, 0.0));
166	0	const int horizontal_x2 = static_cast<int>(
167	0	std::min<double>(dhorizontal_x2, nRasterXSize));
168
169		// Fill the horizontal segment (separately from the rest).
170	0	if (bAvoidBurningSamePoints)
171	0	{
172	0	polyInts2[ints2++] = horizontal_x1;
173	0	polyInts2[ints2++] = horizontal_x2;
174	0	}
175	0	else
176	0	{
177	0	pfnScanlineFunc(
178	0	pCBData, y, horizontal_x1, horizontal_x2 - 1,
179	0	dfVariant == nullptr ? 0 : dfVariant[0]);
180	0	}
181	0	}
182		// else: Skip top horizontal segments.
183		// They are already filled in the regular loop.
184	0	continue;
185	0	}
186
187	0	if (dy < dy2 && dy >= dy1)
188	0	{
189	0	const double intersect = std::clamp<double>(
190	0	(dy - dy1) * (dx2 - dx1) / (dy2 - dy1) + dx1, INT_MIN,
191	0	INT_MAX);
192
193	0	polyInts[ints++] = static_cast<int>(floor(intersect + 0.5));
194	0	}
195	0	}
196
197	0	std::sort(polyInts.begin(), polyInts.begin() + ints);
198	0	std::sort(polyInts2.begin(), polyInts2.begin() + ints2);
199
200	0	for (int i = 0; i + 1 < ints; i += 2)
201	0	{
202	0	if (polyInts[i] <= maxx && polyInts[i + 1] > minx)
203	0	{
204	0	pfnScanlineFunc(pCBData, y, polyInts[i], polyInts[i + 1] - 1,
205	0	dfVariant == nullptr ? 0 : dfVariant[0]);
206	0	}
207	0	}
208
209	0	for (int i2 = 0, i = 0; i2 + 1 < ints2; i2 += 2)
210	0	{
211	0	if (polyInts2[i2] <= maxx && polyInts2[i2 + 1] > minx)
212	0	{
213		// "synchronize" polyInts[i] with polyInts2[i2]
214	0	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	0	if (i + 1 >= ints \|\| polyInts[i] != polyInts2[i2])
219	0	{
220	0	pfnScanlineFunc(pCBData, y, polyInts2[i2],
221	0	polyInts2[i2 + 1] - 1,
222	0	dfVariant == nullptr ? 0 : dfVariant[0]);
223	0	}
224	0	}
225	0	}
226	0	}
227	0	}
228
229		/************************************************************************/
230		/* GDALdllImagePoint() */
231		/************************************************************************/
232
233		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	0	{
238	0	for (int i = 0; i < nPartCount; i++)
239	0	{
240	0	const double dfX = padfX[i];
241	0	const double dfY = padfY[i];
242	0	double dfVariant = 0.0;
243	0	if (padfVariant != nullptr)
244	0	dfVariant = padfVariant[i];
245
246	0	if (0 <= dfX && dfX < nRasterXSize && 0 <= dfY && dfY < nRasterYSize)
247	0	pfnPointFunc(pCBData, static_cast<int>(dfY), static_cast<int>(dfX),
248	0	dfVariant);
249	0	}
250	0	}
251
252		/************************************************************************/
253		/* GDALdllImageLine() */
254		/************************************************************************/
255
256		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	0	{
261	0	if (!nPartCount)
262	0	return;
263
264	0	for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
265	0	{
266	0	for (int j = 1; j < panPartSize[i]; j++)
267	0	{
268	0	double dfX = padfX[n + j - 1];
269	0	double dfY = padfY[n + j - 1];
270	0	double dfXEnd = padfX[n + j];
271	0	double dfYEnd = padfY[n + j];
272
273		// Skip segments that are off the target region.
274	0	if ((dfY < 0.0 && dfYEnd < 0.0) \|\|
275	0	(dfY > nRasterYSize && dfYEnd > nRasterYSize) \|\|
276	0	(dfX < 0.0 && dfXEnd < 0.0) \|\|
277	0	(dfX > nRasterXSize && dfXEnd > nRasterXSize))
278	0	continue;
279
280		// TODO: clamp coordinates to [0, nRasterXSize] * [0, nRasterYSize]
281	0	if (!(dfX >= INT_MIN && dfX <= INT_MAX && dfY >= INT_MIN &&
282	0	dfY <= INT_MAX && dfXEnd >= INT_MIN && dfXEnd <= INT_MAX &&
283	0	dfYEnd >= INT_MIN && dfYEnd <= INT_MAX))
284	0	{
285	0	CPLErrorOnce(
286	0	CE_Warning, CPLE_AppDefined,
287	0	"GDALdllImageLine(): coordinates outside of int range");
288	0	continue;
289	0	}
290
291	0	int iX = static_cast<int>(floor(dfX));
292	0	int iY = static_cast<int>(floor(dfY));
293
294	0	const int iX1 = static_cast<int>(floor(dfXEnd));
295	0	const int iY1 = static_cast<int>(floor(dfYEnd));
296
297	0	double dfVariant = 0.0;
298	0	double dfVariant1 = 0.0;
299	0	if (padfVariant != nullptr &&
300	0	pCBData->eBurnValueSource != GBV_UserBurnValue)
301	0	{
302	0	dfVariant = padfVariant[n + j - 1];
303	0	dfVariant1 = padfVariant[n + j];
304	0	}
305
306	0	int nDeltaX = std::abs(iX1 - iX);
307	0	int nDeltaY = std::abs(iY1 - iY);
308
309		// Step direction depends on line direction.
310	0	const int nXStep = (iX > iX1) ? -1 : 1;
311	0	const int nYStep = (iY > iY1) ? -1 : 1;
312
313		// Determine the line slope.
314	0	if (nDeltaX >= nDeltaY)
315	0	{
316	0	const int nXError = nDeltaY << 1;
317	0	const int nYError = nXError - (nDeltaX << 1);
318	0	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	0	const double dfDeltaVariant =
323	0	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	0	if (j != panPartSize[i] - 1)
329	0	{
330	0	nDeltaX--;
331	0	}
332
333	0	while (nDeltaX-- >= 0)
334	0	{
335	0	if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
336	0	iY < nRasterYSize)
337	0	pfnPointFunc(pCBData, iY, iX, dfVariant);
338
339	0	dfVariant += dfDeltaVariant;
340	0	iX += nXStep;
341	0	if (nError > 0)
342	0	{
343	0	iY += nYStep;
344	0	nError += nYError;
345	0	}
346	0	else
347	0	{
348	0	nError += nXError;
349	0	}
350	0	}
351	0	}
352	0	else
353	0	{
354	0	const int nXError = nDeltaX << 1;
355	0	const int nYError = nXError - (nDeltaY << 1);
356	0	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	0	double dfDeltaVariant =
361	0	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	0	if (j != panPartSize[i] - 1)
367	0	{
368	0	nDeltaY--;
369	0	}
370
371	0	while (nDeltaY-- >= 0)
372	0	{
373	0	if (0 <= iX && iX < nRasterXSize && 0 <= iY &&
374	0	iY < nRasterYSize)
375	0	pfnPointFunc(pCBData, iY, iX, dfVariant);
376
377	0	dfVariant += dfDeltaVariant;
378	0	iY += nYStep;
379	0	if (nError > 0)
380	0	{
381	0	iX += nXStep;
382	0	nError += nYError;
383	0	}
384	0	else
385	0	{
386	0	nError += nXError;
387	0	}
388	0	}
389	0	}
390	0	}
391	0	}
392	0	}
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		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	0	{
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	0	constexpr double EPSILON_INTERSECT_ONLY = 1e-4;
419
420	0	if (!nPartCount)
421	0	return;
422
423	0	for (int i = 0, n = 0; i < nPartCount; n += panPartSize[i++])
424	0	{
425	0	std::set<std::pair<int, int>> lastBurntPoints;
426	0	std::set<std::pair<int, int>> newBurntPoints;
427
428	0	for (int j = 1; j < panPartSize[i]; j++)
429	0	{
430	0	lastBurntPoints = std::move(newBurntPoints);
431	0	newBurntPoints.clear();
432
433	0	double dfX = padfX[n + j - 1];
434	0	double dfY = padfY[n + j - 1];
435
436	0	double dfXEnd = padfX[n + j];
437	0	double dfYEnd = padfY[n + j];
438
439		// Skip segments that are off the target region.
440	0	if ((dfY < 0.0 && dfYEnd < 0.0) \|\|
441	0	(dfY > nRasterYSize && dfYEnd > nRasterYSize) \|\|
442	0	(dfX < 0.0 && dfXEnd < 0.0) \|\|
443	0	(dfX > nRasterXSize && dfXEnd > nRasterXSize))
444	0	continue;
445
446		// TODO: clamp coordinates to [0, nRasterXSize] * [0, nRasterYSize]
447	0	if (!(dfX >= INT_MIN && dfX <= INT_MAX && dfY >= INT_MIN &&
448	0	dfY <= INT_MAX && dfXEnd >= INT_MIN && dfXEnd <= INT_MAX &&
449	0	dfYEnd >= INT_MIN && dfYEnd <= INT_MAX))
450	0	{
451	0	CPLErrorOnce(CE_Warning, CPLE_AppDefined,
452	0	"GDALdllImageLineAllTouched(): coordinates "
453	0	"outside of int range");
454	0	continue;
455	0	}
456
457	0	double dfVariant = 0.0;
458	0	double dfVariantEnd = 0.0;
459	0	if (padfVariant != nullptr &&
460	0	static_cast<GDALRasterizeInfo *>(pCBData)->eBurnValueSource !=
461	0	GBV_UserBurnValue)
462	0	{
463	0	dfVariant = padfVariant[n + j - 1];
464	0	dfVariantEnd = padfVariant[n + j];
465	0	}
466
467		// Swap if needed so we can proceed from left2right (X increasing)
468	0	if (dfX > dfXEnd)
469	0	{
470	0	std::swap(dfX, dfXEnd);
471	0	std::swap(dfY, dfYEnd);
472	0	std::swap(dfVariant, dfVariantEnd);
473	0	}
474
475		// Special case for vertical lines.
476
477	0	if (fabs(dfX - dfXEnd) < .01)
478	0	{
479	0	if (bIntersectOnly)
480	0	{
481	0	if (std::abs(dfX - std::round(dfX)) <
482	0	EPSILON_INTERSECT_ONLY &&
483	0	std::abs(dfXEnd - std::round(dfXEnd)) <
484	0	EPSILON_INTERSECT_ONLY)
485	0	continue;
486	0	}
487
488	0	if (dfYEnd < dfY)
489	0	{
490	0	std::swap(dfY, dfYEnd);
491	0	std::swap(dfVariant, dfVariantEnd);
492	0	}
493
494	0	const int iX = static_cast<int>(floor(dfXEnd));
495	0	int iY = static_cast<int>(floor(dfY));
496	0	int iYEnd =
497	0	static_cast<int>(floor(dfYEnd - EPSILON_INTERSECT_ONLY));
498
499	0	if (iX < 0 \|\| iX >= nRasterXSize)
500	0	continue;
501
502	0	double dfDeltaVariant = 0.0;
503	0	if (dfYEnd - dfY > 0.0)
504	0	dfDeltaVariant = (dfVariantEnd - dfVariant) /
505	0	(dfYEnd - dfY); // Per unit change in iY.
506
507		// Clip to the borders of the target region.
508	0	if (iY < 0)
509	0	iY = 0;
510	0	if (iYEnd >= nRasterYSize)
511	0	iYEnd = nRasterYSize - 1;
512	0	dfVariant += dfDeltaVariant * (iY - dfY);
513
514	0	if (padfVariant == nullptr)
515	0	{
516	0	for (; iY <= iYEnd; iY++)
517	0	{
518	0	if (bAvoidBurningSamePoints)
519	0	{
520	0	auto yx = std::pair<int, int>(iY, iX);
521	0	if (lastBurntPoints.find(yx) !=
522	0	lastBurntPoints.end())
523	0	{
524	0	continue;
525	0	}
526	0	newBurntPoints.insert(yx);
527	0	}
528	0	pfnPointFunc(pCBData, iY, iX, 0.0);
529	0	}
530	0	}
531	0	else
532	0	{
533	0	for (; iY <= iYEnd; iY++, dfVariant += dfDeltaVariant)
534	0	{
535	0	if (bAvoidBurningSamePoints)
536	0	{
537	0	auto yx = std::pair<int, int>(iY, iX);
538	0	if (lastBurntPoints.find(yx) !=
539	0	lastBurntPoints.end())
540	0	{
541	0	continue;
542	0	}
543	0	newBurntPoints.insert(yx);
544	0	}
545	0	pfnPointFunc(pCBData, iY, iX, dfVariant);
546	0	}
547	0	}
548
549	0	continue; // Next segment.
550	0	}
551
552	0	const double dfDeltaVariant =
553	0	(dfVariantEnd - dfVariant) /
554	0	(dfXEnd - dfX); // Per unit change in iX.
555
556		// Special case for horizontal lines.
557	0	if (fabs(dfY - dfYEnd) < .01)
558	0	{
559	0	if (bIntersectOnly)
560	0	{
561	0	if (std::abs(dfY - std::round(dfY)) <
562	0	EPSILON_INTERSECT_ONLY &&
563	0	std::abs(dfYEnd - std::round(dfYEnd)) <
564	0	EPSILON_INTERSECT_ONLY)
565	0	continue;
566	0	}
567
568	0	if (dfXEnd < dfX)
569	0	{
570	0	std::swap(dfX, dfXEnd);
571	0	std::swap(dfVariant, dfVariantEnd);
572	0	}
573
574	0	int iX = static_cast<int>(floor(dfX));
575	0	const int iY = static_cast<int>(floor(dfY));
576	0	int iXEnd =
577	0	static_cast<int>(floor(dfXEnd - EPSILON_INTERSECT_ONLY));
578
579	0	if (iY < 0 \|\| iY >= nRasterYSize)
580	0	continue;
581
582		// Clip to the borders of the target region.
583	0	if (iX < 0)
584	0	iX = 0;
585	0	if (iXEnd >= nRasterXSize)
586	0	iXEnd = nRasterXSize - 1;
587	0	dfVariant += dfDeltaVariant * (iX - dfX);
588
589	0	if (padfVariant == nullptr)
590	0	{
591	0	for (; iX <= iXEnd; iX++)
592	0	{
593	0	if (bAvoidBurningSamePoints)
594	0	{
595	0	auto yx = std::pair<int, int>(iY, iX);
596	0	if (lastBurntPoints.find(yx) !=
597	0	lastBurntPoints.end())
598	0	{
599	0	continue;
600	0	}
601	0	newBurntPoints.insert(yx);
602	0	}
603	0	pfnPointFunc(pCBData, iY, iX, 0.0);
604	0	}
605	0	}
606	0	else
607	0	{
608	0	for (; iX <= iXEnd; iX++, dfVariant += dfDeltaVariant)
609	0	{
610	0	if (bAvoidBurningSamePoints)
611	0	{
612	0	auto yx = std::pair<int, int>(iY, iX);
613	0	if (lastBurntPoints.find(yx) !=
614	0	lastBurntPoints.end())
615	0	{
616	0	continue;
617	0	}
618	0	newBurntPoints.insert(yx);
619	0	}
620	0	pfnPointFunc(pCBData, iY, iX, dfVariant);
621	0	}
622	0	}
623
624	0	continue; // Next segment.
625	0	}
626
627		/* --------------------------------------------------------------------
628		*/
629		/* General case - left to right sloped. */
630		/* --------------------------------------------------------------------
631		*/
632	0	const double dfSlope = (dfYEnd - dfY) / (dfXEnd - dfX);
633
634		// Clip segment in X.
635	0	if (dfXEnd > nRasterXSize)
636	0	{
637	0	dfYEnd -= (dfXEnd - nRasterXSize) * dfSlope;
638	0	dfXEnd = nRasterXSize;
639	0	}
640	0	if (dfX < 0.0)
641	0	{
642	0	dfY += (0.0 - dfX) * dfSlope;
643	0	dfVariant += dfDeltaVariant * (0.0 - dfX);
644	0	dfX = 0.0;
645	0	}
646
647		// Clip segment in Y.
648	0	if (dfYEnd > dfY)
649	0	{
650	0	if (dfY < 0.0)
651	0	{
652	0	const double dfDiffX = (0.0 - dfY) / dfSlope;
653	0	dfX += dfDiffX;
654	0	dfVariant += dfDeltaVariant * dfDiffX;
655	0	dfY = 0.0;
656	0	}
657	0	if (dfYEnd >= nRasterYSize)
658	0	{
659	0	dfXEnd += (dfYEnd - nRasterYSize) / dfSlope;
660	0	if (dfXEnd > nRasterXSize)
661	0	dfXEnd = nRasterXSize;
662		// dfYEnd is no longer used afterwards, but for
663		// consistency it should be:
664		// dfYEnd = nRasterXSize;
665	0	}
666	0	}
667	0	else
668	0	{
669	0	if (dfY >= nRasterYSize)
670	0	{
671	0	const double dfDiffX = (nRasterYSize - dfY) / dfSlope;
672	0	dfX += dfDiffX;
673	0	dfVariant += dfDeltaVariant * dfDiffX;
674	0	dfY = nRasterYSize;
675	0	}
676	0	if (dfYEnd < 0.0)
677	0	{
678	0	dfXEnd -= (dfYEnd - 0) / dfSlope;
679		// dfYEnd is no longer used afterwards, but for
680		// consistency it should be:
681		// dfYEnd = 0.0;
682	0	}
683	0	}
684
685		// Step from pixel to pixel.
686	0	while (dfX >= 0.0 && dfX < dfXEnd)
687	0	{
688	0	const int iX = static_cast<int>(floor(dfX));
689	0	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	0	if (iY >= 0 && iY < nRasterYSize)
695	0	{
696	0	if (bAvoidBurningSamePoints)
697	0	{
698	0	auto yx = std::pair<int, int>(iY, iX);
699	0	if (lastBurntPoints.find(yx) == lastBurntPoints.end() &&
700	0	newBurntPoints.find(yx) == newBurntPoints.end())
701	0	{
702	0	newBurntPoints.insert(yx);
703	0	pfnPointFunc(pCBData, iY, iX, dfVariant);
704	0	}
705	0	}
706	0	else
707	0	{
708	0	pfnPointFunc(pCBData, iY, iX, dfVariant);
709	0	}
710	0	}
711
712	0	double dfStepX = floor(dfX + 1.0) - dfX;
713	0	double dfStepY = dfStepX * dfSlope;
714
715		// Step to right pixel without changing scanline?
716	0	if (static_cast<int>(floor(dfY + dfStepY)) == iY)
717	0	{
718	0	dfX += dfStepX;
719	0	dfY += dfStepY;
720	0	dfVariant += dfDeltaVariant * dfStepX;
721	0	}
722	0	else if (dfSlope < 0)
723	0	{
724	0	dfStepY = iY - dfY;
725	0	if (dfStepY > -0.000000001)
726	0	dfStepY = -0.000000001;
727
728	0	dfStepX = dfStepY / dfSlope;
729	0	dfX += dfStepX;
730	0	dfY += dfStepY;
731	0	dfVariant += dfDeltaVariant * dfStepX;
732	0	}
733	0	else
734	0	{
735	0	dfStepY = (iY + 1) - dfY;
736	0	if (dfStepY < 0.000000001)
737	0	dfStepY = 0.000000001;
738
739	0	dfStepX = dfStepY / dfSlope;
740	0	dfX += dfStepX;
741	0	dfY += dfStepY;
742	0	dfVariant += dfDeltaVariant * dfStepX;
743	0	}
744	0	} // Next step along segment.
745	0	} // Next segment.
746	0	} // Next part.
747	0	}