/* -----------------------------------------------------------------------------

The copyright in this software is being made available under the Clear BSD

License, included below. No patent rights, trademark rights and/or 

other Intellectual Property Rights other than the copyrights concerning 

the Software are granted under this license.

The Clear BSD License

Copyright (c) 2019-2026, Fraunhofer-Gesellschaft zur Förderung der angewandten Forschung e.V. & The VVenC Authors.

All rights reserved.

Redistribution and use in source and binary forms, with or without modification,

are permitted (subject to the limitations in the disclaimer below) provided that

the following conditions are met:

     * Redistributions of source code must retain the above copyright notice,

     this list of conditions and the following disclaimer.

     * Redistributions in binary form must reproduce the above copyright

     notice, this list of conditions and the following disclaimer in the

     documentation and/or other materials provided with the distribution.

     * Neither the name of the copyright holder nor the names of its

     contributors may be used to endorse or promote products derived from this

     software without specific prior written permission.

NO EXPRESS OR IMPLIED LICENSES TO ANY PARTY'S PATENT RIGHTS ARE GRANTED BY

THIS LICENSE. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND

CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A

PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR

CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,

EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,

PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR

BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER

IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)

ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE

POSSIBILITY OF SUCH DAMAGE.

------------------------------------------------------------------------------------------- */

/** \file     MCTF.cpp

\brief    MCTF class

*/

#include "MCTF.h"

#include <math.h>

#include "CommonLib/Picture.h"

#include "CommonLib/dtrace_buffer.h"

#include "Utilities/NoMallocThreadPool.h"

namespace vvenc {

#ifdef TRACE_ENABLE_ITT

static __itt_string_handle* itt_handle_est = __itt_string_handle_create( "MCTF_est" );

static __itt_domain* itt_domain_MCTF_est   = __itt_domain_create( "MCTFEst" );

static __itt_string_handle* itt_handle_flt = __itt_string_handle_create( "MCTF_flt" );

static __itt_domain* itt_domain_MCTF_flt   = __itt_domain_create( "MCTFFlt" );

#endif

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

// Constructor / destructor / initialization / destroy

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

const double MCTF::m_chromaFactor     =  0.55;

const double MCTF::m_sigmaMultiplier  =  9.0;

const int MCTF::m_range               = VVENC_MCTF_RANGE-2;

const int MCTF::m_motionVectorFactor  = 16;

const int MCTF::m_padding             = MCTF_PADDING;

const int16_t MCTF::m_interpolationFilter8[16][8] =

{

  {   0,   0,   0,  64,   0,   0,   0,   0 },   //0

  {   0,   1,  -3,  64,   4,  -2,   0,   0 },   //1 -->-->

  {   0,   1,  -6,  62,   9,  -3,   1,   0 },   //2 -->

  {   0,   2,  -8,  60,  14,  -5,   1,   0 },   //3 -->-->

  {   0,   2,  -9,  57,  19,  -7,   2,   0 },   //4

  {   0,   3, -10,  53,  24,  -8,   2,   0 },   //5 -->-->

  {   0,   3, -11,  50,  29,  -9,   2,   0 },   //6 -->

  {   0,   3, -11,  44,  35, -10,   3,   0 },   //7 -->-->

  {   0,   1,  -7,  38,  38,  -7,   1,   0 },   //8

  {   0,   3, -10,  35,  44, -11,   3,   0 },   //9 -->-->

  {   0,   2,  -9,  29,  50, -11,   3,   0 },   //10-->

  {   0,   2,  -8,  24,  53, -10,   3,   0 },   //11-->-->

  {   0,   2,  -7,  19,  57,  -9,   2,   0 },   //12

  {   0,   1,  -5,  14,  60,  -8,   2,   0 },   //13-->-->

  {   0,   1,  -3,   9,  62,  -6,   1,   0 },   //14-->

  {   0,   0,  -2,   4,  64,  -3,   1,   0 }    //15-->-->

};

const int16_t MCTF::m_interpolationFilter4[16][4] =

{

  {  0, 64,  0,  0 },    //0

  { -2, 62,  4,  0 },    //1 -->-->

  { -2, 58, 10, -2 },    //2 -->

  { -4, 56, 14, -2 },    //3 -->-->

  { -4, 54, 16, -2 },    //4

  { -6, 52, 20, -2 },    //5 -->-->

  { -6, 46, 28, -4 },    //6 -->

  { -4, 42, 30, -4 },    //7 -->-->

  { -4, 36, 36, -4 },    //8

  { -4, 30, 42, -4 },    //9 -->-->

  { -4, 28, 46, -6 },    //10-->

  { -2, 20, 52, -6 },    //11-->-->

  { -2, 16, 54, -4 },    //12

  { -2, 14, 56, -4 },    //13-->-->

  { -2, 10, 58, -2 },    //14-->

  {  0,  4, 62, -2 },    //15-->-->

};

const double MCTF::m_refStrengths[2][6] = // min(..., {3 or 5} / (1 + 2 * |POC offset|))

{ // abs(POC offset)

  // 1       2       3       4       5       6

  { 0.84375, 0.6, 0.4286, 0.3333, 0.2727, 0.2308 }, // RA

  { 1.12500, 1.0, 0.7143, 0.5556, 0.4545, 0.3846 }  // LD

};

const int    MCTF::m_cuTreeThresh[4] = { 75, 60,     30, 15 };

const double MCTF::m_cuTreeCenter    =           45;

int motionErrorLumaInt( const Pel* org, const ptrdiff_t origStride, const Pel* buf, const ptrdiff_t buffStride, const int w, const int h, const int besterror )

{

  int error = 0;

  for( int y1 = 0; y1 < h; y1++ )

  {

    const Pel* origRowStart   = org + y1 * origStride;

    const Pel* bufferRowStart = buf + y1 * buffStride;

    for( int x1 = 0; x1 < w; x1 += 2 )

    {

      int diff = origRowStart[x1] - bufferRowStart[x1];

      error += diff * diff;

      diff = origRowStart[x1 + 1] - bufferRowStart[x1 + 1];

      error += diff * diff;

    }

    if( error > besterror )

    {

      return error;

    }

  }

  return error;

}

int motionErrorLumaFrac6( const Pel *org, const ptrdiff_t origStride, const Pel *buf, const ptrdiff_t buffStride, const int w, const int h, const int16_t *xFilter, const int16_t *yFilter, const int bitDepth, const int besterror )

{

  int error = 0;

  Pel tempArray[64 + 8][64];

  int sum, base;

  const Pel maxSampleValue = ( 1 << bitDepth ) - 1;

  for( int y1 = 1; y1 < h + 7; y1++ )

  {

    const int yOffset = y1 - 3;

    const Pel *sourceRow = buf + yOffset * buffStride;

    for( int x1 = 0; x1 < w; x1++ )

    {

      sum = 0;

      base = x1 - 3;

      const Pel *rowStart = sourceRow + base;

      sum += xFilter[1] * rowStart[1];

      sum += xFilter[2] * rowStart[2];

      sum += xFilter[3] * rowStart[3];

      sum += xFilter[4] * rowStart[4];

      sum += xFilter[5] * rowStart[5];

      sum += xFilter[6] * rowStart[6];

      sum = ( sum + ( 1 << 5 ) ) >> 6;

      sum = sum < 0 ? 0 : ( sum > maxSampleValue ? maxSampleValue : sum );

      tempArray[y1][x1] = sum;

    }

  }

  for( int y1 = 0; y1 < h; y1++ )

  {

    const Pel *origRow = org + y1 * origStride;

    for( int x1 = 0; x1 < w; x1++ )

    {

      sum = 0;

      sum += yFilter[1] * tempArray[y1 + 1][x1];

      sum += yFilter[2] * tempArray[y1 + 2][x1];

      sum += yFilter[3] * tempArray[y1 + 3][x1];

      sum += yFilter[4] * tempArray[y1 + 4][x1];

      sum += yFilter[5] * tempArray[y1 + 5][x1];

      sum += yFilter[6] * tempArray[y1 + 6][x1];

      sum = ( sum + ( 1 << 5 ) ) >> 6;

      sum = sum < 0 ? 0 : ( sum > maxSampleValue ? maxSampleValue : sum );

      error += ( sum - origRow[x1] ) * ( sum - origRow[x1] );

    }

    if( error > besterror )

    {

      return error;

    }

  }

  return error;

}

int motionErrorLumaFrac4( const Pel* org, const ptrdiff_t origStride, const Pel* buf, const ptrdiff_t buffStride, const int w, const int h, const int16_t* xFilter, const int16_t* yFilter, const int bitDepth, const int besterror )

{

  int error = 0;

  Pel tempArray[64 + 4][64];

  int sum, base;

  const Pel maxSampleValue = ( 1 << bitDepth ) - 1;

  for( int y1 = 0; y1 < h + 3; y1++ )

  {

    const int yOffset = y1 - 1;

    const Pel* sourceRow = buf + yOffset * buffStride;

    for( int x1 = 0; x1 < w; x1++ )

    {

      sum = 0;

      base = x1 - 1;

      const Pel* rowStart = sourceRow + base;

      sum += xFilter[0] * rowStart[0];

      sum += xFilter[1] * rowStart[1];

      sum += xFilter[2] * rowStart[2];

      sum += xFilter[3] * rowStart[3];

      sum = ( sum + ( 1 << 5 ) ) >> 6;

      sum = sum < 0 ? 0 : ( sum > maxSampleValue ? maxSampleValue : sum );

      tempArray[y1][x1] = sum;

    }

  }

  for( int y1 = 0; y1 < h; y1++ )

  {

    const Pel* origRow = org + y1 * origStride;

    for( int x1 = 0; x1 < w; x1++ )

    {

      sum = 0;

      sum += yFilter[0] * tempArray[y1 + 0][x1];

      sum += yFilter[1] * tempArray[y1 + 1][x1];

      sum += yFilter[2] * tempArray[y1 + 2][x1];

      sum += yFilter[3] * tempArray[y1 + 3][x1];

      sum = ( sum + ( 1 << 5 ) ) >> 6;

      sum = sum < 0 ? 0 : ( sum > maxSampleValue ? maxSampleValue : sum );

      error += ( sum - origRow[x1] ) * ( sum - origRow[x1] );

    }

    if( error > besterror )

    {

      return error;

    }

  }

  return error;

}

void applyFrac8Core_6Tap( const Pel* org, const ptrdiff_t origStride, Pel* dst, const ptrdiff_t dstStride, const int w, const int h, const int16_t* xFilter, const int16_t* yFilter, const int bitDepth )

{

  const int numFilterTaps   = 7;

  const int centreTapOffset = 3;

  const int maxValue        = ( 1 << bitDepth ) - 1;

  Pel tempArray[64 + numFilterTaps][64];

  for( int by = 1; by < h + numFilterTaps - 1; by++ )

  {

    const int yOffset = by - centreTapOffset;

    const Pel *sourceRow = org + yOffset * origStride;

    for( int bx = 0; bx < w; bx++ )

    {

      int base = bx - centreTapOffset;

      const Pel *rowStart = sourceRow + base;

      int sum = 0;

      sum += xFilter[1] * rowStart[1];

      sum += xFilter[2] * rowStart[2];

      sum += xFilter[3] * rowStart[3];

      sum += xFilter[4] * rowStart[4];

      sum += xFilter[5] * rowStart[5];

      sum += xFilter[6] * rowStart[6];

      sum = ( sum + ( 1 << 5 ) ) >> 6;

      tempArray[by][bx] = sum;

    }

  }

  Pel *dstRow = dst;

  for( int by = 0; by < h; by++, dstRow += dstStride )

  {

    Pel *dstPel = dstRow;

    for( int bx = 0; bx < w; bx++, dstPel++ )

    {

      int sum = 0;

      sum += yFilter[1] * tempArray[by + 1][bx];

      sum += yFilter[2] * tempArray[by + 2][bx];

      sum += yFilter[3] * tempArray[by + 3][bx];

      sum += yFilter[4] * tempArray[by + 4][bx];

      sum += yFilter[5] * tempArray[by + 5][bx];

      sum += yFilter[6] * tempArray[by + 6][bx];

      sum = ( sum + ( 1 << 5 ) ) >> 6;

      sum = sum < 0 ? 0 : ( sum > maxValue ? maxValue : sum );

      *dstPel = sum;

    }

  }

}

void applyFrac8Core_4Tap( const Pel* org, const ptrdiff_t origStride, Pel* dst, const ptrdiff_t dstStride, const int w, const int h, const int16_t* xFilter, const int16_t* yFilter, const int bitDepth )

{

  const int numFilterTaps   = 3;

  const int centreTapOffset = 1;

  const int maxValue        = ( 1 << bitDepth ) - 1;

  Pel tempArray[64 + numFilterTaps][64];

  for( int by = 0; by < h + numFilterTaps; by++ )

  {

    const int yOffset    = by - centreTapOffset;

    const Pel* sourceRow = org + yOffset * origStride;

    for( int bx = 0; bx < w; bx++ )

    {

      int base = bx - centreTapOffset;

      const Pel* rowStart = sourceRow + base;

      int sum = 0;

      sum += xFilter[0] * rowStart[0];

      sum += xFilter[1] * rowStart[1];

      sum += xFilter[2] * rowStart[2];

      sum += xFilter[3] * rowStart[3];

      sum = ( sum + ( 1 << 5 ) ) >> 6;

      tempArray[by][bx] = sum;

    }

  }

  Pel* dstRow = dst;

  for( int by = 0; by < h; by++, dstRow += dstStride )

  {

    Pel* dstPel = dstRow;

    for( int bx = 0; bx < w; bx++, dstPel++ )

    {

      int sum = 0;

      sum += yFilter[0] * tempArray[by + 0][bx];

      sum += yFilter[1] * tempArray[by + 1][bx];

      sum += yFilter[2] * tempArray[by + 2][bx];

      sum += yFilter[3] * tempArray[by + 3][bx];

      sum = ( sum + ( 1 << 5 ) ) >> 6;

      sum = sum < 0 ? 0 : ( sum > maxValue ? maxValue : sum );

      *dstPel = sum;

    }

  }

}

inline static float fastExp( float n, float d )

{

  // using the e^x ~= ( 1 + x/n )^n for n -> inf

  float x = 1.0f + n / ( d * 1024 );

  x *= x; x *= x; x *= x; x *= x;

  x *= x; x *= x; x *= x; x *= x;

  x *= x; x *= x;

  return x;

}

static const int32_t xSzm[6] = {0, 1, 20, 336, 5440, 87296};

// works for bit depths up to incl. 12 and power-of-2 block dimensions in both directions

void applyPlanarCorrectionCore( const Pel* refPel, const ptrdiff_t refStride, Pel* dstPel, const ptrdiff_t dstStride, const int32_t w, const int32_t h, const ClpRng& clpRng, const uint16_t motionError )

{

  const int32_t blockSize = w * h;

  const int32_t log2Width = floorLog2 (w);

  const int32_t maxPelVal = clpRng.max();

  const int32_t mWeight   = std::min (512u, (uint32_t) motionError * (uint32_t) motionError);

  const int32_t xSum      = (blockSize * (w - 1)) >> 1;

  int32_t x1yzm = 0,  x2yzm = 0,  ySum = 0;

  int32_t b0, b1, b2;

  int64_t numer, denom;

  for (int32_t y = 0; y < h; y++) // sum up dot-products between indices and sample diffs

  {

    for (int32_t x = 0; x < w; x++)

    {

      const Pel* pDst = dstPel + y * dstStride + x;

      const Pel* pRef = refPel + y * refStride + x;

      const int32_t z = *pDst - *pRef;

      x1yzm += x * z;  x2yzm += y * z;  ySum += z;

    }

  }

  denom = blockSize * xSzm[log2Width]; // plane-fit parameters, in fixed-point arithmetic

  numer = (int64_t) mWeight * ((int64_t) x1yzm * blockSize - xSum * ySum);

  b1 = int32_t ((numer < 0 ? numer - (denom >> 1) : numer + (denom >> 1)) / denom);

  b1 = (b1 < INT16_MIN ? INT16_MIN : (b1 > INT16_MAX ? INT16_MAX : b1));

  numer = (int64_t) mWeight * ((int64_t) x2yzm * blockSize - xSum * ySum);

  b2 = int32_t ((numer < 0 ? numer - (denom >> 1) : numer + (denom >> 1)) / denom);

  b2 = (b2 > INT16_MAX ? INT16_MAX : (b2 < INT16_MIN ? INT16_MIN : b2));

  b0 = (mWeight * ySum - (b1 + b2) * xSum + (blockSize >> 1)) >> (log2Width << 1);

  if (b0 == 0 && b1 == 0 && b2 == 0) return;

  for (int32_t y = 0; y < h; y++) // perform deblocking by adding fitted correction plane

  {

    for (int32_t x = 0; x < w; x++)

    {

      Pel* const pDst = dstPel + y * dstStride + x;

      const int32_t p = (b0 + b1 * x + b2 * y + 256) >> 9; // fixed-point plane corrector

      const int32_t z = *pDst - p;

      *pDst = Pel (z < 0 ? 0 : (z > maxPelVal ? maxPelVal : z));

    }

  }

}

void applyBlockCore( const CPelBuf& src, PelBuf& dst, const CompArea& blk, const ClpRng& clpRng, const Pel** correctedPics, int numRefs, const int* verror, const double* refStrenghts, double weightScaling, double sigmaSq )

{

  const int         w = blk.width;

  const int         h = blk.height;

  const int        bx = blk.x;

  const int        by = blk.y;

  const ptrdiff_t srcStride = src.stride;

  const ptrdiff_t dstStride = dst.stride;

  const Pel *srcPel = src.bufAt( bx, by );

        Pel *dstPel = dst.bufAt( bx, by );

  const Pel maxSampleValue = clpRng.max();

  int vnoise[2 * VVENC_MCTF_RANGE] = { 0, };

  float vsw [2 * VVENC_MCTF_RANGE] = { 0.0f, };

  float vww [2 * VVENC_MCTF_RANGE] = { 0.0f, };

  int minError = INT32_MAX;

  for( int i = 0; i < numRefs; i++ )

  {

    int64_t variance = 0, diffsum = 0;

    const ptrdiff_t refStride = w;

    const Pel *     refPel    = correctedPics[i];

    for( int y1 = 0; y1 < h; y1++ )

    {

      for( int x1 = 0; x1 < w; x1++ )

      {

        const Pel pix = *( srcPel + srcStride * y1 + x1 );

        const Pel ref = *( refPel + refStride * y1 + x1 );

        const int diff = pix - ref;

        variance += diff * diff;

        if( x1 != w - 1 )

        {

          const Pel pixR = *( srcPel + srcStride * y1 + x1 + 1 );

          const Pel refR = *( refPel + refStride * y1 + x1 + 1 );

          const int diffR = pixR - refR;

          diffsum += ( diffR - diff ) * ( diffR - diff );

        }

        if( y1 != h - 1 )

        {

          const Pel pixD = *( srcPel + srcStride * y1 + x1 + srcStride );

          const Pel refD = *( refPel + refStride * y1 + x1 + refStride );

          const int diffD = pixD - refD;

          diffsum += ( diffD - diff ) * ( diffD - diff );

        }

      }

    }

    variance *= (int64_t) 1 << (2*(10-clpRng.bd));

    diffsum  *= (int64_t) 1 << (2*(10-clpRng.bd));

    const int cntV = w * h;

    const int cntD = 2 * cntV - w - h;

    vnoise[i] = ( int ) round( ( 15.0 * cntD / cntV * variance + 5.0 ) / ( diffsum + 5.0 ) );

    minError = std::min( minError, verror[i] );

  }

  for( int i = 0; i < numRefs; i++ )

  {

    const int error = verror[i];

    const int noise = vnoise[i];

    float ww = 1, sw = 1;

    ww *= ( noise < 25 ) ? 1.0 : 0.6;

    sw *= ( noise < 25 ) ? 1.0 : 0.8;

    ww *= ( error < 50 ) ? 1.2 : ( ( error > 100 ) ? 0.6 : 1.0 );

    sw *= ( error < 50 ) ? 1.0 : 0.8;

    ww *= ( ( minError + 1.0 ) / ( error + 1.0 ) );

    vww[i] = ww * weightScaling * refStrenghts[i];

    vsw[i] = sw * 2 * sigmaSq;

  }

  for( int y = 0; y < h; y++ )

  {

    for( int x = 0; x < w; x++ )

    {

      const Pel orgVal  = *( srcPel + srcStride * y + x );

      float temporalWeightSum = 1.0;

      float newVal = ( float ) orgVal;

      for( int i = 0; i < numRefs; i++ )

      {

        const Pel* pCorrectedPelPtr = correctedPics[i] + y * w + x;

        const int    refVal = *pCorrectedPelPtr;

        const int    diff   = refVal - orgVal;

        const float  diffSq = diff * diff;

        float weight = vww[i] * fastExp( -diffSq, vsw[i] );

        newVal += weight * refVal;

        temporalWeightSum += weight;

      }

      newVal /= temporalWeightSum;

      Pel sampleVal = ( Pel ) ( newVal + 0.5 );

      sampleVal = ( sampleVal < 0 ? 0 : ( sampleVal > maxSampleValue ? maxSampleValue : sampleVal ) );

      *( dstPel + dstStride * y + x ) = sampleVal;

    }

  }

}

double calcVarCore( const Pel* org, const ptrdiff_t origStride, const int w, const int h )

{

  // calculate average

  int avg = 0;

  for( int y1 = 0; y1 < h; y1++ )

  {

    for( int x1 = 0; x1 < w; x1++ )

    {

      avg = avg + *( org + x1 + y1 * origStride );

    }

  }

  avg <<= 4;

  avg = avg / ( w * h );

  // calculate variance

  int64_t variance = 0;

  for( int y1 = 0; y1 < h; y1++ )

  {

    for( int x1 = 0; x1 < w; x1++ )

    {

      int pix = *( org + x1 + y1 * origStride ) << 4;

      variance = variance + ( pix - avg ) * ( pix - avg );

    }

  }

  return variance / 256.0;

}

MCTF::MCTF( bool enableOpt )

  : m_encCfg     ( nullptr )

  , m_threadPool ( nullptr )

  , m_isFinalPass( true )

  , m_filterPoc  ( 0 )

  , m_lastPicIn  ( nullptr )

{

  m_motionErrorLumaIntX     = motionErrorLumaInt;

  m_motionErrorLumaInt8     = motionErrorLumaInt;

  m_motionErrorLumaFracX[0] = motionErrorLumaFrac6;

  m_motionErrorLumaFrac8[0] = motionErrorLumaFrac6;

  m_motionErrorLumaFracX[1] = motionErrorLumaFrac4;

  m_motionErrorLumaFrac8[1] = motionErrorLumaFrac4;

  m_applyFrac[0][0]         = applyFrac8Core_6Tap;

  m_applyFrac[0][1]         = applyFrac8Core_4Tap;

  m_applyFrac[1][0]         = applyFrac8Core_6Tap;

  m_applyFrac[1][1]         = applyFrac8Core_4Tap;

  m_applyPlanarCorrection   = applyPlanarCorrectionCore;

  m_applyBlock              = applyBlockCore;

  m_calcVar                 = calcVarCore;

  if( enableOpt )

  {

#if defined( TARGET_SIMD_X86 ) && ENABLE_SIMD_OPT_MCTF

    initMCTF_X86();

#endif

#if defined( TARGET_SIMD_ARM ) && ENABLE_SIMD_OPT_MCTF

    initMCTF_ARM();

#endif

  }

}

MCTF::~MCTF()

{

}

void MCTF::init( const VVEncCfg& encCfg, bool isFinalPass, NoMallocThreadPool* threadPool )

{

  CHECK( encCfg.m_vvencMCTF.numFrames != encCfg.m_vvencMCTF.numStrength, "should have been checked before" );

  m_encCfg      = &encCfg;

  m_threadPool  = threadPool;

  m_isFinalPass = isFinalPass;

  m_filterPoc   = 0;

  m_area        = Area( 0, 0, m_encCfg->m_PadSourceWidth, m_encCfg->m_PadSourceHeight );

  // TLayer (TL) dependent definition of drop frames: TL = 4,  TL = 3,  TL = 2,  TL = 1,  TL = 0

  const static int sMCTFSpeed[5] { 0, 0, ((3<<12) + (2<<9) + (2<<6) + (0<<3) + 0),   ((3<<12) + (2<<9) + (2<<6) + (0<<3) + 0),   ((3<<12) + (3<<9) + (3<<6) + (2<<3) + 2) };

  m_MCTFSpeedVal     = sMCTFSpeed[ m_encCfg->m_vvencMCTF.MCTFSpeed ];

  m_lowResFltSearch  = m_encCfg->m_vvencMCTF.MCTFSpeed > 0;

  m_searchPttrn      = m_encCfg->m_vvencMCTF.MCTFSpeed > 0 ? ( m_encCfg->m_vvencMCTF.MCTFSpeed >= 3 ? 2 : 1 ) : 0;

  m_mctfUnitSize     = m_encCfg->m_vvencMCTF.MCTFUnitSize;

}

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

// Public member functions

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

void MCTF::initPicture( Picture* pic )

{

  pic->getOrigBuf().extendBorderPel( MCTF_FRAME_MARGIN, true );

  pic->setSccFlags( m_encCfg );

}

void MCTF::processPictures( const PicList& picList, AccessUnitList& auList, PicList& doneList, PicList& freeList )

{

  // ensure this is only processed if necessary 

  if( picList.empty() || ( m_lastPicIn == picList.back() && ! picList.back()->isFlush ))

  {

    return;

  }

  m_lastPicIn = picList.back();

  // filter one picture (either all or up to frames to be encoded)

  if( picList.size()

      && m_filterPoc <= picList.back()->poc

      && ( m_encCfg->m_framesToBeEncoded <= 0 || m_filterPoc < m_encCfg->m_framesToBeEncoded ) )

  {

    // setup fifo of pictures to be filtered

    std::deque<Picture*> picFifo;

    int filterIdx = 0;

    for( auto pic : picList )

    {

      const int minPoc = m_filterPoc - VVENC_MCTF_RANGE;

      const int maxPoc = m_encCfg->m_vvencMCTF.MCTFFutureReference ? m_filterPoc + VVENC_MCTF_RANGE : m_filterPoc;

      if( pic->poc >= minPoc && pic->poc <= maxPoc )

      {

        picFifo.push_back( pic );

        if( pic->poc < m_filterPoc )

        {

          filterIdx += 1;

        }

      }

    }

    CHECK( picFifo.empty(), "MCTF: no pictures to be filtered found" );

    CHECK( filterIdx >= (int)picFifo.size(), "MCTF: picture filter error" );

    CHECK( picFifo[ filterIdx ]->poc != m_filterPoc, "MCTF: picture filter error" );

    // filter picture (when more than 1 picture is available for processing)

    if( picFifo.size() > 1 )

    {

      filter( picFifo, filterIdx );

    }

    // set picture done

    doneList.push_back( picFifo[ filterIdx ] );

  }

  // mark pictures not needed anymore

  for( auto pic : picList )

  {

    if( pic->poc > m_filterPoc - VVENC_MCTF_RANGE )

      break;

    freeList.push_back( pic );

  }

  m_filterPoc += 1;

}

void MCTF::motionEstimationMCTF(Picture* curPic, std::deque<TemporalFilterSourcePicInfo> &srcFrameInfo, const PelStorage& origBuf, PelStorage& origSubsampled2, PelStorage& origSubsampled4, PelStorage& origSubsampled8, std::vector<double> &mvErr, double &minError, bool addLevel, bool calcErr)

{

  srcFrameInfo.push_back(TemporalFilterSourcePicInfo());

  TemporalFilterSourcePicInfo& srcPic = srcFrameInfo.back();

  const int wInBlks = (m_area.width + m_mctfUnitSize - 1) / m_mctfUnitSize;

  const int hInBlks = (m_area.height + m_mctfUnitSize - 1) / m_mctfUnitSize;

  srcPic.picBuffer.createFromBuf(curPic->getOrigBuf());

  srcPic.mvs.allocate(wInBlks, hInBlks);

  srcPic.index = std::min(5, std::abs(curPic->poc - m_filterPoc) - 1);

  {

    const int width = m_area.width;

    const int height = m_area.height;

    Array2D<MotionVector> mv_0(width / (m_mctfUnitSize * 8) + 1, height / (m_mctfUnitSize * 8) + 1);

    Array2D<MotionVector> mv_1(width / (m_mctfUnitSize * 4) + 1, height / (m_mctfUnitSize * 4) + 1);

    Array2D<MotionVector> mv_2(width / (m_mctfUnitSize * 2) + 1, height / (m_mctfUnitSize * 2) + 1);

    PelStorage bufferSub2;

    PelStorage bufferSub4;

    subsampleLuma(srcPic.picBuffer, bufferSub2);

    subsampleLuma(bufferSub2, bufferSub4);

    if (addLevel)

    {

      Array2D<MotionVector> mv_m(width / (m_mctfUnitSize * 16) + 1, height / (m_mctfUnitSize * 16) + 1);

      PelStorage bufferSub8;

      subsampleLuma(bufferSub4, bufferSub8);

      motionEstimationLuma(mv_m, origSubsampled8, bufferSub8, m_padding, 2 * m_mctfUnitSize);

      motionEstimationLuma(mv_0, origSubsampled4, bufferSub4, m_padding, 2 * m_mctfUnitSize, &mv_m, 2);

    }

    else

    {

      motionEstimationLuma(mv_0, origSubsampled4, bufferSub4, m_padding, 2 * m_mctfUnitSize);

    }

    motionEstimationLuma(mv_1, origSubsampled2, bufferSub2, m_padding, 2 * m_mctfUnitSize, &mv_0, 2);

    motionEstimationLuma(mv_2, origBuf, srcPic.picBuffer, MCTF_FRAME_MARGIN, 2 * m_mctfUnitSize, &mv_1, 2);

    motionEstimationLuma(srcPic.mvs, origBuf, srcPic.picBuffer, MCTF_FRAME_MARGIN, m_mctfUnitSize, &mv_2, 1, true);

    if (calcErr)

    {

      double sumErr = 0.0;

      for (int y = 0; y < srcPic.mvs.h(); y++) // going over ref pic in block steps

      {

        for (int x = 0; x < srcPic.mvs.w(); x++)

        {

          sumErr += srcPic.mvs.get(x, y).error;

        }

      }

      double S = 1.0 / (srcPic.mvs.w() * srcPic.mvs.h());

      mvErr.push_back(sumErr * S);

      minError = std::min(minError, sumErr * S);

    }

  }

}

void MCTF::filter( const std::deque<Picture*>& picFifo, int filterIdx )

{

  PROFILER_SCOPE_AND_STAGE( 1, g_timeProfiler, P_MCTF );

  Picture* pic = picFifo[ filterIdx ];

  // first-pass temporal downsampling

  if( ! m_isFinalPass && pic->gopEntry->m_skipFirstPass )

  {

    return;

  }

  const int mctfIdx            = pic->gopEntry ? pic->gopEntry->m_mctfIndex : -1;

  const double overallStrength = mctfIdx >= 0 ? m_encCfg->m_vvencMCTF.MCTFStrengths[ mctfIdx ] : -1.0;

  double   meanRmsAcrossPic    = 0.0;

  uint64_t sumSRmsAcrossPic    = 0;

  uint16_t nMax = 0, maxRmsCTU = 0;

  bool  isFilterThisFrame      = mctfIdx >= 0;

  int dropFrames = ( m_encCfg->m_usePerceptQPA ? VVENC_MCTF_RANGE >> 1 : 0 );

  if( mctfIdx >= 0 )

  {

    const int idxTLayer = m_encCfg->m_vvencMCTF.numFrames - (mctfIdx + 1);

    const int threshold = (m_MCTFSpeedVal >> (idxTLayer * 3)) & 7;

    dropFrames          = std::min(VVENC_MCTF_RANGE, threshold);

    isFilterThisFrame   = threshold < VVENC_MCTF_RANGE;

  }

  const int filterFrames = VVENC_MCTF_RANGE - 2 - dropFrames;

  int dropFramesFront = std::min( std::max(                                          filterIdx - filterFrames, 0 ), dropFrames + 2 );

  int dropFramesBack  = std::min( std::max( static_cast<int>( picFifo.size() ) - 1 - filterIdx - filterFrames, 0 ), dropFrames + 2 );

  if( !pic->useMCTF && !pic->gopEntry->m_isStartOfGop )

  {

    isFilterThisFrame = false;

  }

  if ( isFilterThisFrame )

  {

    bool  useMCTFadaptation = true;

    const bool condAddLevel = useMCTFadaptation && m_area.width >= 1920;

    std::vector<double> mvErr;

    double minError = MAX_DOUBLE;

    const PelStorage& origBuf = pic->getOrigBuffer();

          PelStorage& fltrBuf = pic->getFilteredOrigBuffer();

    // subsample original picture so it only needs to be done once

    PelStorage origSubsampled2;

    PelStorage origSubsampled4;

    PelStorage origSubsampled8;

    subsampleLuma( origBuf,         origSubsampled2 );

    subsampleLuma( origSubsampled2, origSubsampled4 );

    if (condAddLevel)

    {

      subsampleLuma(origSubsampled4, origSubsampled8);

    }

    // determine motion vectors

    std::deque<TemporalFilterSourcePicInfo> srcFrameInfo;

    for ( int i = dropFramesFront; i < picFifo.size() - dropFramesBack; i++ )

    {

      Picture* curPic = picFifo[ i ];

      if ( curPic->poc == m_filterPoc )

      {

        continue;

      }

      motionEstimationMCTF(curPic, srcFrameInfo, origBuf, origSubsampled2, origSubsampled4, origSubsampled8 ,mvErr, minError, condAddLevel, useMCTFadaptation);

    }

    int lastIndexRefFr = -1;

    if ((m_encCfg->m_vvencMCTF.MCTFSpeed < 4) && (minError > 80))

    {

      useMCTFadaptation = false;

    }

    if (useMCTFadaptation && minError)

    {

      const double errThr = 0.75 * minError * srcFrameInfo.size();

      int avgErrCond = 0;

      int minErrCond = 0;

      double factErr = m_encCfg->m_vvencMCTF.MCTFSpeed < 4 ? 1.0 : 2.0 ;

      double SizeThi = m_encCfg->m_vvencMCTF.MCTFSpeed < 4  ? filterFrames + 1 : 3.0;

      for (const double& framMvErr : mvErr)

      {

        if (factErr * framMvErr > errThr)

        {

          avgErrCond++;

        }

        if (framMvErr > SizeThi * minError)

        {

          minErrCond++;

        }

      }

      int newFilterFrames = minErrCond ? filterFrames : (filterFrames + 2 - avgErrCond);

      if (filterFrames <= 2 && newFilterFrames > 3)   newFilterFrames = 3;

      for (int curIdx = filterFrames + 1; (curIdx < newFilterFrames + 1)&&((lastIndexRefFr == -1)); curIdx++)

      {

        for (int i = 0; i < picFifo.size(); i++)

        {

          Picture* curPic = picFifo[i];

          if (curIdx == std::abs(curPic->poc - m_filterPoc))

          {

            motionEstimationMCTF(curPic, srcFrameInfo, origBuf, origSubsampled2, origSubsampled4, origSubsampled8, mvErr, minError, condAddLevel, m_encCfg->m_vvencMCTF.MCTFSpeed == 4);

            if (m_encCfg->m_vvencMCTF.MCTFSpeed == 4)

            {

              int nSize = (int(srcFrameInfo.size()) & 1) + int(srcFrameInfo.size());

              const double errThrcur = 0.75 * minError * nSize;

              if (mvErr.back() > errThrcur)

              {

                lastIndexRefFr = curIdx;

                break;

              }

            }

          }

        }

      }

      if ((lastIndexRefFr != -1))

      {

        for (auto it = srcFrameInfo.begin(); it != srcFrameInfo.end(); )

        {

          if ((it->index + 1) >= lastIndexRefFr)

          {

            it = srcFrameInfo.erase(it);

          }

          else

          {

            ++it;

          }

        }

      }

    }

    // filter

    if( pic->useMCTF )

    {

      fltrBuf.create( m_encCfg->m_internChromaFormat, m_area, 0, MCTF_FRAME_MARGIN );

      bilateralFilter( origBuf, srcFrameInfo, fltrBuf, overallStrength );

    }

    if( m_encCfg->m_blockImportanceMapping || m_encCfg->m_usePerceptQPA || pic->gopEntry->m_isStartOfGop )

    {

      const int ctuSize        = m_encCfg->m_bimCtuSize;

      const int widthInCtus    = ( m_area.width  + ctuSize - 1 ) / ctuSize;

      const int heightInCtus   = ( m_area.height + ctuSize - 1 ) / ctuSize;

      const int numCtu         = widthInCtus * heightInCtus;

      const int ctuBlocks      = ctuSize / m_mctfUnitSize;

      std::vector<double> sumError( numCtu * 2, 0 );

      std::vector<uint32_t> sumRMS( numCtu * 2, 0 ); // RMS of motion estimation error

      std::vector<uint16_t> maxRMS( numCtu * 2, 0 ); // maximum block estimation error

      std::vector<double> blkCount( numCtu * 2, 0 );

      int distFactor[2] = { 3,3 };

      for( auto& srcPic : srcFrameInfo )

      {

        if( srcPic.index >= 2 )

        {