#include "trop.h"

#include "loop_macros.h"

#include "tpixelutils.h"

#include "quickputP.h"

#ifndef TNZCORE_LIGHT

#include "tpalette.h"

#include "tcolorstyles.h"

#endif

/*

#ifndef __sgi

#include <algorithm>

#endif

*/

// The following must be old IRIX code. Should be re-tested.

// It seems that gcc compiles it, but requiring a LOT of

// resources... very suspect...

/*#ifdef __LP64__

#include "optimize_for_lp64.h"

#endif*/

//=============================================================================

//=============================================================================

//=============================================================================

#ifdef OPTIMIZE_FOR_LP64

void quickResample_optimized(const TRasterP &dn, const TRasterP &up,

                             const TAffine &aff,

                             TRop::ResampleFilterType filterType);

#endif

namespace {

inline TPixel32 applyColorScale(const TPixel32 &color,

                                const TPixel32 &colorScale,

                                bool toBePremultiplied = false) {

  /*--

   * Prevent colors from being darkened when quickputting a semi-transparent

   * raster on the Viewer

   * --*/

  if (colorScale.r == 0 && colorScale.g == 0 && colorScale.b == 0) {

    /*--

     * When toBePremultiplied is ON, Premultiply is done later, so it is not

     * done here.

     * --*/

    if (toBePremultiplied)

      return TPixel32(color.r, color.g, color.b, color.m * colorScale.m / 255);

    else

      return TPixel32(

          color.r * colorScale.m / 255, color.g * colorScale.m / 255,

          color.b * colorScale.m / 255, color.m * colorScale.m / 255);

  }

  int r = 255 - (255 - color.r) * (255 - colorScale.r) / 255;

  int g = 255 - (255 - color.g) * (255 - colorScale.g) / 255;

  int b = 255 - (255 - color.b) * (255 - colorScale.b) / 255;

  return premultiply(TPixel32(r, g, b, color.m * colorScale.m / 255));

}

//------------------------------------------------------------------------------

inline TPixel32 applyColorScaleCMapped(const TPixel32 &color,

                                       const TPixel32 &colorScale) {

  int r = 255 - (255 - color.r) * (255 - colorScale.r) / 255;

  int g = 255 - (255 - color.g) * (255 - colorScale.g) / 255;

  int b = 255 - (255 - color.b) * (255 - colorScale.b) / 255;

  return premultiply(TPixel32(r, g, b, color.m * colorScale.m / 255));

}

//------------------------------------------------------------------------------

void doQuickPutFilter(const TRaster32P &dn, const TRaster32P &up,

                      const TAffine &aff) {

  //  se aff e' degenere la controimmagine di up e' un segmento (o un punto)

  if ((aff.a11 * aff.a22 - aff.a12 * aff.a21) == 0) return;

  //  contatore bit di shift

  const int PADN = 16;

  //  maschera del filtro bilineare

  const int MASKN = (1 << PADN) - 1;

  //  max dimensioni di up gestibili (limite imposto dal numero di bit

  //  disponibili per la parte intera di xL, yL)

  assert(std::max(up->getLx(), up->getLy()) <

         (1 << (8 * sizeof(int) - PADN - 1)));

  TRectD boundingBoxD = TRectD(convert(dn->getSize())) *

                        (aff * TRectD(0, 0, up->getLx() - 2, up->getLy() - 2));

  //  clipping

  if (boundingBoxD.x0 >= boundingBoxD.x1 || boundingBoxD.y0 >= boundingBoxD.y1)

    return;

  //  clipping y su dn

  int yMin = std::max(tfloor(boundingBoxD.y0), 0);

  //  clipping y su dn

  int yMax = std::min(tceil(boundingBoxD.y1), dn->getLy() - 1);

  //  clipping x su dn

  int xMin = std::max(tfloor(boundingBoxD.x0), 0);

  //  clipping x su dn

  int xMax = std::min(tceil(boundingBoxD.x1), dn->getLx() - 1);

  TAffine invAff = inv(aff);  //  inversa di aff

  //  nel disegnare la y-esima scanline di dn, il passaggio al pixel successivo

  //  comporta l'incremento (deltaXD, deltaYD) delle coordinate del pixel

  //  corrispondente di up

  double deltaXD = invAff.a11;

  double deltaYD = invAff.a21;

  //  deltaXD "TLonghizzato" (round)

  int deltaXL = tround(deltaXD * (1 << PADN));

  //  deltaYD "TLonghizzato" (round)

  int deltaYL = tround(deltaYD * (1 << PADN));

  //  se aff "TLonghizzata" (round) e' degenere la controimmagine di up e' un

  // segmento (o un punto)

  if ((deltaXL == 0) && (deltaYL == 0)) return;

  //	naturale predecessore di up->getLx() - 1

  int lxPred = (up->getLx() - 2) * (1 << PADN);

  //	naturale predecessore di up->getLy() - 1

  int lyPred = (up->getLy() - 2) * (1 << PADN);

  int dnWrap = dn->getWrap();

  int upWrap = up->getWrap();

  dn->lock();

  up->lock();

  TPixel32 *dnRow     = dn->pixels(yMin);

  TPixel32 *upBasePix = up->pixels();

  //  scorre le scanline di boundingBoxD

  for (int y = yMin; y <= yMax; y++, dnRow += dnWrap) {

    //	(1)    equazione k-parametrica della y-esima

    //               scanline di boundingBoxD:

    //	       (xMin, y) + k*(1, 0),  k = 0, ..., (xMax - xMin)

    //	(2)    equazione k-parametrica dell'immagine mediante

    //               invAff di (1):

    //	       invAff*(xMin, y) + k*(deltaXD, deltaYD),

    //	       k = kMin, ..., kMax

    //               con 0 <= kMin <= kMax <= (xMax - xMin)

    //  calcola kMin, kMax per la scanline corrente intersecando la (2)

    //  con i lati di up

    //  il segmento [a, b] di up e' la controimmagine mediante aff

    //  della porzione di scanline  [ (xMin, y), (xMax, y) ] di dn

    //        TPointD b = invAff*TPointD(xMax, y);

    TPointD a = invAff * TPointD(xMin, y);

    //	(xL0, yL0) sono le coordinate di a in versione "TLonghizzata"

    //	0 <= xL0 + k*deltaXL

    //          <= (up->getLx() - 2)*(1<

    //          <= kMinX

    //          <= kMin

    //          <= k

    //          <= kMax

    //          <= kMaxX

    //          <= (xMax - xMin)

    //

    //	0 <= yL0 + k*deltaYL

    //          <= (up->getLy() - 2)*(1<

    //          <= kMinY

    //          <= kMin

    //          <= k

    //          <= kMax

    //          <= kMaxY

    //          <= (yMax - yMin)

    int xL0 = tround(a.x * (1 << PADN));  //  xL0 inizializzato

    int yL0 = tround(a.y * (1 << PADN));  //  yL0 inizializzato

    //  calcola kMinX, kMaxX, kMinY, kMaxY

    int kMinX = 0, kMaxX = xMax - xMin;  //  clipping su dn

    int kMinY = 0, kMaxY = yMax - yMin;  //  clipping su dn

    //        0 <= xL0 + k*deltaXL <= (up->getLx() - 2)*(1<<PADN)

    //                   <=>

    //        0 <= xL0 + k*deltaXL <= lxPred

    //

    //

    //	0 <= yL0 + k*deltaYL <= (up->getLy() - 2)*(1<<PADN)

    //                   <=>

    //        0 <= yL0 + k*deltaYL <= lyPred

    //  calcola kMinX, kMaxX

    if (deltaXL == 0) {

      //  [a, b] verticale esterno ad up contratto

      if ((xL0 < 0) || (lxPred < xL0)) continue;

      //  altrimenti usa solo

      //  kMinY, kMaxY ((deltaXL != 0) || (deltaYL != 0))

    } else if (deltaXL > 0) {

      if (lxPred < xL0)  //  [a, b] esterno ad up+(bordo destro)

        continue;

      kMaxX = (lxPred - xL0) / deltaXL;  //  floor

      if (xL0 < 0) {

        kMinX = ((-xL0) + deltaXL - 1) / deltaXL;  //  ceil

      }

    } else  //  (deltaXL < 0)

    {

      if (xL0 < 0)  //  [a, b] esterno ad up contratto

        continue;

      kMaxX = xL0 / (-deltaXL);  //  floor

      if (lxPred < xL0) {

        kMinX = (xL0 - lxPred - deltaXL - 1) / (-deltaXL);  //  ceil

      }

    }

    //  calcola kMinY, kMaxY

    if (deltaYL == 0) {

      //  [a, b] orizzontale esterno ad up contratto

      if ((yL0 < 0) || (lyPred < yL0)) continue;

      //  altrimenti usa solo

      //  kMinX, kMaxX ((deltaXL != 0) || (deltaYL != 0))

    } else if (deltaYL > 0) {

      if (lyPred < yL0)  //  [a, b] esterno ad up contratto

        continue;

      kMaxY = (lyPred - yL0) / deltaYL;  //  floor

      if (yL0 < 0) {

        kMinY = ((-yL0) + deltaYL - 1) / deltaYL;  //  ceil

      }

    } else  //  (deltaYL < 0)

    {

      if (yL0 < 0)  //  [a, b] esterno ad up contratto

        continue;

      kMaxY = yL0 / (-deltaYL);  //  floor

      if (lyPred < yL0) {

        kMinY = (yL0 - lyPred - deltaYL - 1) / (-deltaYL);  //  ceil

      }

    }

    //  calcola kMin, kMax effettuando anche il clipping su dn

    int kMin = std::max({kMinX, kMinY, (int)0});

    int kMax = std::min({kMaxX, kMaxY, xMax - xMin});

    TPixel32 *dnPix    = dnRow + xMin + kMin;

    TPixel32 *dnEndPix = dnRow + xMin + kMax + 1;

    //  (xL, yL) sono le coordinate (inizializzate per il round)

    //  in versione "TLonghizzata"

    //	del pixel corrente di up

    int xL = xL0 + (kMin - 1) * deltaXL;  //  inizializza xL

    int yL = yL0 + (kMin - 1) * deltaYL;  //  inizializza yL

    //  scorre i pixel sulla y-esima scanline di boundingBoxD

    for (; dnPix < dnEndPix; ++dnPix) {

      xL += deltaXL;

      yL += deltaYL;

      //  il punto di up TPointD(xL/(1<

      //  e' approssimato con (xI, yI)

      int xI = xL >> PADN;  //	troncato

      int yI = yL >> PADN;  //	troncato

      assert((0 <= xI) && (xI <= up->getLx() - 1) && (0 <= yI) &&

             (yI <= up->getLy() - 1));

      //  (xI, yI)

      TPixel32 *upPix00 = upBasePix + (yI * upWrap + xI);

      //  (xI + 1, yI)

      TPixel32 *upPix10 = upPix00 + 1;

      //  (xI, yI + 1)

      TPixel32 *upPix01 = upPix00 + upWrap;

      //  (xI + 1, yI + 1)

      TPixel32 *upPix11 = upPix00 + upWrap + 1;

      //  filtro bilineare 4 pixels: calcolo dei pesi

      int xWeight1 = (xL & MASKN);

      int xWeight0 = (1 << PADN) - xWeight1;

      int yWeight1 = (yL & MASKN);

      int yWeight0 = (1 << PADN) - yWeight1;

      //  filtro bilineare 4 pixels: media pesata sui singoli canali

      int rColDownTmp =

          (xWeight0 * (upPix00->r) + xWeight1 * ((upPix10)->r)) >> PADN;

      int gColDownTmp =

          (xWeight0 * (upPix00->g) + xWeight1 * ((upPix10)->g)) >> PADN;

      int bColDownTmp =

          (xWeight0 * (upPix00->b) + xWeight1 * ((upPix10)->b)) >> PADN;

      int rColUpTmp =

          (xWeight0 * ((upPix01)->r) + xWeight1 * ((upPix11)->r)) >> PADN;

      int gColUpTmp =

          (xWeight0 * ((upPix01)->g) + xWeight1 * ((upPix11)->g)) >> PADN;

      int bColUpTmp =

          (xWeight0 * ((upPix01)->b) + xWeight1 * ((upPix11)->b)) >> PADN;

      unsigned char rCol =

          (unsigned char)((yWeight0 * rColDownTmp + yWeight1 * rColUpTmp) >>

                          PADN);

      unsigned char gCol =

          (unsigned char)((yWeight0 * gColDownTmp + yWeight1 * gColUpTmp) >>

                          PADN);

      unsigned char bCol =

          (unsigned char)((yWeight0 * bColDownTmp + yWeight1 * bColUpTmp) >>

                          PADN);

      TPixel32 upPix = TPixel32(rCol, gCol, bCol, upPix00->m);

      if (upPix.m == 0)

        continue;

      else if (upPix.m == 255)

        *dnPix = upPix;

      else

        *dnPix = quickOverPix(*dnPix, upPix);

    }

  }

  dn->unlock();

  up->unlock();

}

//=============================================================================

//=============================================================================

//=============================================================================

void doQuickPutNoFilter(const TRaster32P &dn, const TRaster32P &up,

                        const TAffine &aff, const TPixel32 &colorScale,

                        bool doPremultiply, bool whiteTransp, bool firstColumn,

                        bool doRasterDarkenBlendedView) {

  //  se aff := TAffine(sx, 0, tx, 0, sy, ty) e' degenere la controimmagine

  //  di up e' un segmento (o un punto)

  if ((aff.a11 * aff.a22 - aff.a12 * aff.a21) == 0) return;

  //  contatore bit di shift

  const int PADN = 16;

  //  max dimensioni di up gestibili (limite imposto dal numero di bit

  //  disponibili per la parte intera di xL, yL)

  assert(std::max(up->getLx(), up->getLy()) <

         (1 << (8 * sizeof(int) - PADN - 1)));

  TRectD boundingBoxD =

      TRectD(convert(dn->getBounds())) *

      (aff * TRectD(-0.5, -0.5, up->getLx() - 0.5, up->getLy() - 0.5));

  //  clipping

  if (boundingBoxD.x0 >= boundingBoxD.x1 || boundingBoxD.y0 >= boundingBoxD.y1)

    return;

  //  clipping y su dn

  int yMin = std::max(tfloor(boundingBoxD.y0), 0);

  //  clipping y su dn

  int yMax = std::min(tceil(boundingBoxD.y1), dn->getLy() - 1);

  //  clipping x su dn

  int xMin = std::max(tfloor(boundingBoxD.x0), 0);

  //  clipping x su dn

  int xMax = std::min(tceil(boundingBoxD.x1), dn->getLx() - 1);

  //  inversa di aff

  TAffine invAff = inv(aff);

  //  nel disegnare la y-esima scanline di dn, il passaggio al pixel

  //  successivo comporta l'incremento (deltaXD, deltaYD) delle coordinate del

  //  pixel corrispondente di up

  double deltaXD = invAff.a11;

  double deltaYD = invAff.a21;

  //  deltaXD "TLonghizzato" (round)

  int deltaXL = tround(deltaXD * (1 << PADN));

  //  deltaYD "TLonghizzato" (round)

  int deltaYL = tround(deltaYD * (1 << PADN));

  //  se aff "TLonghizzata" (round) e' degenere la controimmagine di up e' un

  //  segmento (o un punto)

  if ((deltaXL == 0) && (deltaYL == 0)) return;

  //  TINT32 predecessore di up->getLx()

  int lxPred = up->getLx() * (1 << PADN) - 1;

  //  TINT32 predecessore di up->getLy()

  int lyPred = up->getLy() * (1 << PADN) - 1;

  int dnWrap = dn->getWrap();

  int upWrap = up->getWrap();

  dn->lock();

  up->lock();

  TPixel32 *dnRow     = dn->pixels(yMin);

  TPixel32 *upBasePix = up->pixels();

  //  scorre le scanline di boundingBoxD

  for (int y = yMin; y <= yMax; y++, dnRow += dnWrap) {

    //  (1)  equazione k-parametrica della y-esima scanline di boundingBoxD:

    //       (xMin, y) + k*(1, 0),  k = 0, ..., (xMax - xMin)

    //  (2)  equazione k-parametrica dell'immagine mediante invAff di (1):

    //       invAff*(xMin, y) + k*(deltaXD, deltaYD),

    //       k = kMin, ..., kMax con 0 <= kMin <= kMax <= (xMax - xMin)

    //  calcola kMin, kMax per la scanline corrente

    //  intersecando la (2) con i lati di up

    //  il segmento [a, b] di up e' la controimmagine mediante aff della

    //  porzione di scanline  [ (xMin, y), (xMax, y) ] di dn

    //  TPointD b = invAff*TPointD(xMax, y);

    TPointD a = invAff * TPointD(xMin, y);

    //  (xL0, yL0) sono le coordinate di a (inizializzate per il round)

    //  in versione "TLonghizzata"

    //  0 <= xL0 + k*deltaXL

    //    <  up->getLx()*(1<<PADN)

    //

    //  0 <= kMinX

    //    <= kMin

    //    <= k

    //    <= kMax

    //    <= kMaxX

    //    <= (xMax - xMin)

    //

    //  0 <= yL0 + k*deltaYL

    //    < up->getLy()*(1<<PADN)

    //  0 <= kMinY

    //    <= kMin

    //    <= k

    //    <= kMax

    //    <= kMaxY

    //    <= (xMax - xMin)

    //  xL0 inizializzato per il round

    int xL0 = tround((a.x + 0.5) * (1 << PADN));

    //  yL0 inizializzato per il round

    int yL0 = tround((a.y + 0.5) * (1 << PADN));

    //  calcola kMinX, kMaxX, kMinY, kMaxY

    int kMinX = 0, kMaxX = xMax - xMin;  //  clipping su dn

    int kMinY = 0, kMaxY = xMax - xMin;  //  clipping su dn

    //  0 <= xL0 + k*deltaXL

    //    < up->getLx()*(1<<PADN)

    //           <=>

    //  0 <= xL0 + k*deltaXL

    //    <= lxPred

    //

    //  0 <= yL0 + k*deltaYL

    //    < up->getLy()*(1<<PADN)

    //           <=>

    //  0 <= yL0 + k*deltaYL

    //    <= lyPred

    //  calcola kMinX, kMaxX

    if (deltaXL == 0) {

      // [a, b] verticale esterno ad up+(bordo destro/basso)

      if ((xL0 < 0) || (lxPred < xL0)) continue;

      //  altrimenti usa solo

      //  kMinY, kMaxY ((deltaXL != 0) || (deltaYL != 0))

    } else if (deltaXL > 0) {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (lxPred < xL0) continue;

      kMaxX = (lxPred - xL0) / deltaXL;  //  floor

      if (xL0 < 0) {

        kMinX = ((-xL0) + deltaXL - 1) / deltaXL;  //  ceil

      }

    } else  //  (deltaXL < 0)

    {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (xL0 < 0) continue;

      kMaxX = xL0 / (-deltaXL);  //  floor

      if (lxPred < xL0) {

        kMinX = (xL0 - lxPred - deltaXL - 1) / (-deltaXL);  //  ceil

      }

    }

    //  calcola kMinY, kMaxY

    if (deltaYL == 0) {

      //  [a, b] orizzontale esterno ad up+(bordo destro/basso)

      if ((yL0 < 0) || (lyPred < yL0)) continue;

      // altrimenti usa solo

      // kMinX, kMaxX ((deltaXL != 0) || (deltaYL != 0))

    } else if (deltaYL > 0) {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (lyPred < yL0) continue;

      kMaxY = (lyPred - yL0) / deltaYL;  //  floor

      if (yL0 < 0) {

        kMinY = ((-yL0) + deltaYL - 1) / deltaYL;  //  ceil

      }

    } else  //  (deltaYL < 0)

    {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (yL0 < 0) continue;

      kMaxY = yL0 / (-deltaYL);  //  floor

      if (lyPred < yL0) {

        kMinY = (yL0 - lyPred - deltaYL - 1) / (-deltaYL);  //  ceil

      }

    }

    //  calcola kMin, kMax effettuando anche il clippind su dn

    int kMin = std::max({kMinX, kMinY, (int)0});

    int kMax = std::min({kMaxX, kMaxY, xMax - xMin});

    TPixel32 *dnPix    = dnRow + xMin + kMin;

    TPixel32 *dnEndPix = dnRow + xMin + kMax + 1;

    //  (xL, yL) sono le coordinate (inizializzate per il round)

    //  in versione "TLonghizzata" del pixel corrente di up

    int xL = xL0 + (kMin - 1) * deltaXL;  //  inizializza xL

    int yL = yL0 + (kMin - 1) * deltaYL;  //  inizializza yL

    //  scorre i pixel sulla y-esima scanline di boundingBoxD

    for (; dnPix < dnEndPix; ++dnPix) {

      xL += deltaXL;

      yL += deltaYL;

      //  il punto di up TPointD(xL/(1<

      //  approssimato con (xI, yI)

      int xI = xL >> PADN;  //  round

      int yI = yL >> PADN;  //  round

      assert((0 <= xI) && (xI <= up->getLx() - 1) && (0 <= yI) &&

             (yI <= up->getLy() - 1));

      TPixel32 upPix = *(upBasePix + (yI * upWrap + xI));

      if (firstColumn) upPix.m = 255;

      if (upPix.m == 0 || (whiteTransp && upPix == TPixel::White)) continue;

      if (colorScale != TPixel32::Black)

        upPix = applyColorScale(upPix, colorScale, doPremultiply);

      if (doRasterDarkenBlendedView)

        *dnPix = quickOverPixDarkenBlended(*dnPix, upPix);

      else {

        if (upPix.m == 255)

          *dnPix = upPix;

        else if (doPremultiply)

          *dnPix = quickOverPixPremult(*dnPix, upPix);

        else

          *dnPix = quickOverPix(*dnPix, upPix);

      }

    }

  }

  dn->unlock();

  up->unlock();

}

//=============================================================================

//=============================================================================

//=============================================================================

void doQuickPutNoFilter(const TRaster32P &dn, const TRaster64P &up,

                        const TAffine &aff, bool doPremultiply,

                        bool firstColumn) {

  //  se aff := TAffine(sx, 0, tx, 0, sy, ty) e' degenere la controimmagine

  //  di up e' un segmento (o un punto)

  if ((aff.a11 * aff.a22 - aff.a12 * aff.a21) == 0) return;

  //  contatore bit di shift

  const int PADN = 16;

  //  max dimensioni di up gestibili (limite imposto dal numero di bit

  //  disponibili per la parte intera di xL, yL)

  assert(std::max(up->getLx(), up->getLy()) <

         (1 << (8 * sizeof(int) - PADN - 1)));

  TRectD boundingBoxD =

      TRectD(convert(dn->getBounds())) *

      (aff * TRectD(-0.5, -0.5, up->getLx() - 0.5, up->getLy() - 0.5));

  //  clipping

  if (boundingBoxD.x0 >= boundingBoxD.x1 || boundingBoxD.y0 >= boundingBoxD.y1)

    return;

  //  clipping y su dn

  int yMin = std::max(tfloor(boundingBoxD.y0), 0);

  //  clipping y su dn

  int yMax = std::min(tceil(boundingBoxD.y1), dn->getLy() - 1);

  //  clipping x su dn

  int xMin = std::max(tfloor(boundingBoxD.x0), 0);

  //  clipping x su dn

  int xMax = std::min(tceil(boundingBoxD.x1), dn->getLx() - 1);

  //  inversa di aff

  TAffine invAff = inv(aff);

  //  nel disegnare la y-esima scanline di dn, il passaggio al pixel

  //  successivo comporta l'incremento (deltaXD, deltaYD) delle coordinate del

  //  pixel corrispondente di up

  double deltaXD = invAff.a11;

  double deltaYD = invAff.a21;

  //  deltaXD "TLonghizzato" (round)

  int deltaXL = tround(deltaXD * (1 << PADN));

  //  deltaYD "TLonghizzato" (round)

  int deltaYL = tround(deltaYD * (1 << PADN));

  //  se aff "TLonghizzata" (round) e' degenere la controimmagine di up e' un

  //  segmento (o un punto)

  if ((deltaXL == 0) && (deltaYL == 0)) return;

  //  TINT32 predecessore di up->getLx()

  int lxPred = up->getLx() * (1 << PADN) - 1;

  //  TINT32 predecessore di up->getLy()

  int lyPred = up->getLy() * (1 << PADN) - 1;

  int dnWrap = dn->getWrap();

  int upWrap = up->getWrap();

  dn->lock();

  up->lock();

  TPixel32 *dnRow     = dn->pixels(yMin);

  TPixel64 *upBasePix = up->pixels();

  //  scorre le scanline di boundingBoxD

  for (int y = yMin; y <= yMax; y++, dnRow += dnWrap) {

    //  (1)  equazione k-parametrica della y-esima scanline di boundingBoxD:

    //       (xMin, y) + k*(1, 0),  k = 0, ..., (xMax - xMin)

    //  (2)  equazione k-parametrica dell'immagine mediante invAff di (1):

    //       invAff*(xMin, y) + k*(deltaXD, deltaYD),

    //       k = kMin, ..., kMax con 0 <= kMin <= kMax <= (xMax - xMin)

    //  calcola kMin, kMax per la scanline corrente

    //  intersecando la (2) con i lati di up

    //  il segmento [a, b] di up e' la controimmagine mediante aff della

    //  porzione di scanline  [ (xMin, y), (xMax, y) ] di dn

    //  TPointD b = invAff*TPointD(xMax, y);

    TPointD a = invAff * TPointD(xMin, y);

    //  (xL0, yL0) sono le coordinate di a (inizializzate per il round)

    //  in versione "TLonghizzata"

    //  0 <= xL0 + k*deltaXL

    //    <  up->getLx()*(1<<PADN)

    //

    //  0 <= kMinX

    //    <= kMin

    //    <= k

    //    <= kMax

    //    <= kMaxX

    //    <= (xMax - xMin)

    //

    //  0 <= yL0 + k*deltaYL

    //    < up->getLy()*(1<<PADN)

    //  0 <= kMinY

    //    <= kMin

    //    <= k

    //    <= kMax

    //    <= kMaxY

    //    <= (xMax - xMin)

    //  xL0 inizializzato per il round

    int xL0 = tround((a.x + 0.5) * (1 << PADN));

    //  yL0 inizializzato per il round

    int yL0 = tround((a.y + 0.5) * (1 << PADN));

    //  calcola kMinX, kMaxX, kMinY, kMaxY

    int kMinX = 0, kMaxX = xMax - xMin;  //  clipping su dn

    int kMinY = 0, kMaxY = xMax - xMin;  //  clipping su dn

    //  0 <= xL0 + k*deltaXL

    //    < up->getLx()*(1<<PADN)

    //           <=>

    //  0 <= xL0 + k*deltaXL

    //    <= lxPred

    //

    //  0 <= yL0 + k*deltaYL

    //    < up->getLy()*(1<<PADN)

    //           <=>

    //  0 <= yL0 + k*deltaYL

    //    <= lyPred

    //  calcola kMinX, kMaxX

    if (deltaXL == 0) {

      // [a, b] verticale esterno ad up+(bordo destro/basso)

      if ((xL0 < 0) || (lxPred < xL0)) continue;

      //  altrimenti usa solo

      //  kMinY, kMaxY ((deltaXL != 0) || (deltaYL != 0))

    } else if (deltaXL > 0) {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (lxPred < xL0) continue;

      kMaxX = (lxPred - xL0) / deltaXL;  //  floor

      if (xL0 < 0) {

        kMinX = ((-xL0) + deltaXL - 1) / deltaXL;  //  ceil

      }

    } else  //  (deltaXL < 0)

    {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (xL0 < 0) continue;

      kMaxX = xL0 / (-deltaXL);  //  floor

      if (lxPred < xL0) {

        kMinX = (xL0 - lxPred - deltaXL - 1) / (-deltaXL);  //  ceil

      }

    }

    //  calcola kMinY, kMaxY

    if (deltaYL == 0) {

      //  [a, b] orizzontale esterno ad up+(bordo destro/basso)

      if ((yL0 < 0) || (lyPred < yL0)) continue;

      // altrimenti usa solo

      // kMinX, kMaxX ((deltaXL != 0) || (deltaYL != 0))

    } else if (deltaYL > 0) {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (lyPred < yL0) continue;

      kMaxY = (lyPred - yL0) / deltaYL;  //  floor

      if (yL0 < 0) {

        kMinY = ((-yL0) + deltaYL - 1) / deltaYL;  //  ceil

      }

    } else  //  (deltaYL < 0)

    {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (yL0 < 0) continue;

      kMaxY = yL0 / (-deltaYL);  //  floor

      if (lyPred < yL0) {

        kMinY = (yL0 - lyPred - deltaYL - 1) / (-deltaYL);  //  ceil

      }

    }

    //  calcola kMin, kMax effettuando anche il clippind su dn

    int kMin = std::max({kMinX, kMinY, (int)0});

    int kMax = std::min({kMaxX, kMaxY, xMax - xMin});

    TPixel32 *dnPix    = dnRow + xMin + kMin;

    TPixel32 *dnEndPix = dnRow + xMin + kMax + 1;

    //  (xL, yL) sono le coordinate (inizializzate per il round)

    //  in versione "TLonghizzata" del pixel corrente di up

    int xL = xL0 + (kMin - 1) * deltaXL;  //  inizializza xL

    int yL = yL0 + (kMin - 1) * deltaYL;  //  inizializza yL

    //  scorre i pixel sulla y-esima scanline di boundingBoxD

    for (; dnPix < dnEndPix; ++dnPix) {

      xL += deltaXL;

      yL += deltaYL;

      //  il punto di up TPointD(xL/(1<

      //  approssimato con (xI, yI)

      int xI = xL >> PADN;  //  round

      int yI = yL >> PADN;  //  round

      assert((0 <= xI) && (xI <= up->getLx() - 1) && (0 <= yI) &&

             (yI <= up->getLy() - 1));

      TPixel64 *upPix = upBasePix + (yI * upWrap + xI);

      if (firstColumn) upPix->m = 65535;

      if (upPix->m == 0)

        continue;

      else if (upPix->m == 65535)

        *dnPix = PixelConverter<TPixel32>::from(*upPix);

      else if (doPremultiply)

        *dnPix =

            quickOverPixPremult(*dnPix, PixelConverter<TPixel32>::from(*upPix));

      else

        *dnPix = quickOverPix(*dnPix, PixelConverter<TPixel32>::from(*upPix));

    }

  }

  dn->unlock();

  up->unlock();

}

//=============================================================================

//=============================================================================

//=============================================================================

void doQuickPutNoFilter(const TRaster32P &dn, const TRasterGR8P &up,

                        const TAffine &aff, const TPixel32 &colorScale) {

  if ((aff.a11 * aff.a22 - aff.a12 * aff.a21) == 0) return;

  const int PADN = 16;

  assert(std::max(up->getLx(), up->getLy()) <

         (1 << (8 * sizeof(int) - PADN - 1)));

  TRectD boundingBoxD =

      TRectD(convert(dn->getBounds())) *

      (aff * TRectD(-0.5, -0.5, up->getLx() - 0.5, up->getLy() - 0.5));

  if (boundingBoxD.x0 >= boundingBoxD.x1 || boundingBoxD.y0 >= boundingBoxD.y1)

    return;

  int yMin = std::max(tfloor(boundingBoxD.y0), 0);

  int yMax = std::min(tceil(boundingBoxD.y1), dn->getLy() - 1);

  int xMin = std::max(tfloor(boundingBoxD.x0), 0);

  int xMax = std::min(tceil(boundingBoxD.x1), dn->getLx() - 1);

  TAffine invAff = inv(aff);

  double deltaXD = invAff.a11;

  double deltaYD = invAff.a21;

  int deltaXL = tround(deltaXD * (1 << PADN));

  int deltaYL = tround(deltaYD * (1 << PADN));

  if ((deltaXL == 0) && (deltaYL == 0)) return;

  int lxPred = up->getLx() * (1 << PADN) - 1;

  int lyPred = up->getLy() * (1 << PADN) - 1;

  int dnWrap = dn->getWrap();

  int upWrap = up->getWrap();

  dn->lock();

  up->lock();

  TPixel32 *dnRow      = dn->pixels(yMin);

  TPixelGR8 *upBasePix = up->pixels();

  for (int y = yMin; y <= yMax; y++, dnRow += dnWrap) {

    TPointD a = invAff * TPointD(xMin, y);

    int xL0 = tround((a.x + 0.5) * (1 << PADN));

    int yL0 = tround((a.y + 0.5) * (1 << PADN));

    int kMinX = 0, kMaxX = xMax - xMin;  //  clipping su dn

    int kMinY = 0, kMaxY = xMax - xMin;  //  clipping su dn

    if (deltaXL == 0) {

      if ((xL0 < 0) || (lxPred < xL0)) continue;

    } else if (deltaXL > 0) {

      if (lxPred < xL0) continue;

      kMaxX = (lxPred - xL0) / deltaXL;  //  floor

      if (xL0 < 0) {

        kMinX = ((-xL0) + deltaXL - 1) / deltaXL;  //  ceil

      }

    } else  //  (deltaXL < 0)

    {

      //  [a, b] esterno ad up+(bordo destro/basso)

      if (xL0 < 0) continue;

      kMaxX = xL0 / (-deltaXL);  //  floor

      if (lxPred < xL0) {

        kMinX = (xL0 - lxPred - deltaXL - 1) / (-deltaXL);  //  ceil

      }

    }

    //  calcola kMinY, kMaxY

    if (deltaYL == 0) {

      //  [a, b] orizzontale esterno ad up+(bordo destro/basso)

      if ((yL0 < 0) || (lyPred < yL0)) continue;

      // altrimenti usa solo

      // kMinX, kMaxX ((deltaXL != 0) || (deltaYL != 0))