Implement StatsLoop

src/ImageSharp/Formats/WebP/BitWriter/Vp8BitWriter.cs

@@ -28,8 +28,6 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
 
         private int maxPos;
 
-        // private bool error;
-
         /// <summary>
         /// Initializes a new instance of the <see cref="Vp8BitWriter"/> class.
         /// </summary>
@@ -43,8 +41,6 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
             this.nbBits = -8;
             this.pos = 0;
             this.maxPos = 0;
-
-            // this.error = false;
         }
 
         /// <inheritdoc/>
@@ -69,13 +65,13 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
                 int v = sign ? -c : c;
                 if (!this.PutBit(v != 0, p.Probabilities[1]))
                 {
-                    p = residual.Prob[WebPConstants.Bands[n]].Probabilities[0];
+                    p = residual.Prob[WebPConstants.Vp8EncBands[n]].Probabilities[0];
                     continue;
                 }
 
                 if (!this.PutBit(v > 1, p.Probabilities[2]))
                 {
-                    p = residual.Prob[WebPConstants.Bands[n]].Probabilities[1];
+                    p = residual.Prob[WebPConstants.Vp8EncBands[n]].Probabilities[1];
                 }
                 else
                 {
@@ -102,7 +98,6 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
                     {
                         int mask;
                         byte[] tab;
-                        var tabIdx = 0;
                         if (v < 3 + (8 << 1))
                         {
                             // VP8Cat3  (3b)
@@ -111,7 +106,6 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
                             v -= 3 + (8 << 0);
                             mask = 1 << 2;
                             tab = WebPConstants.Cat3;
-                            tabIdx = 0;
                         }
                         else if (v < 3 + (8 << 2))
                         {
@@ -121,7 +115,6 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
                             v -= 3 + (8 << 1);
                             mask = 1 << 3;
                             tab = WebPConstants.Cat4;
-                            tabIdx = 0;
                         }
                         else if (v < 3 + (8 << 3))
                         {
@@ -131,7 +124,6 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
                             v -= 3 + (8 << 2);
                             mask = 1 << 4;
                             tab = WebPConstants.Cat5;
-                            tabIdx = 0;
                         }
                         else
                         {
@@ -141,9 +133,9 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
                             v -= 3 + (8 << 3);
                             mask = 1 << 10;
                             tab = WebPConstants.Cat6;
-                            tabIdx = 0;
                         }
 
+                        var tabIdx = 0;
                         while (mask != 0)
                         {
                             this.PutBit(v & mask, tab[tabIdx++]);
@@ -151,7 +143,7 @@ namespace SixLabors.ImageSharp.Formats.WebP.BitWriter
                         }
                     }
 
-                    p = residual.Prob[WebPConstants.Bands[n]].Probabilities[2];
+                    p = residual.Prob[WebPConstants.Vp8EncBands[n]].Probabilities[2];
                 }
 
                 this.PutBitUniform(sign ? 1 : 0);

src/ImageSharp/Formats/WebP/Lossy/PassStats.cs

@@ -0,0 +1,76 @@
+// Copyright (c) Six Labors.
+// Licensed under the Apache License, Version 2.0.
+
+namespace SixLabors.ImageSharp.Formats.WebP.Lossy
+{
+    /// <summary>
+    /// Class for organizing convergence in either size or PSNR.
+    /// </summary>
+    internal class PassStats
+    {
+        public PassStats(long targetSize, float targetPsnr, int qMin, int qMax, int quality)
+        {
+            bool doSizeSearch = targetSize != 0;
+
+            this.IsFirst = true;
+            this.Dq = 10.0f;
+            this.Qmin = qMin;
+            this.Qmax = qMax;
+            this.Q = quality.Clamp(qMin, qMax);
+            this.LastQ = this.Q;
+            this.Target = doSizeSearch ? targetSize
+                : (targetPsnr > 0.0f) ? targetPsnr
+                : 40.0f;   // default, just in case
+            this.Value = 0.0f;
+            this.LastValue = 0.0f;
+            this.DoSizeSearch = doSizeSearch;
+        }
+
+        public bool IsFirst { get; set; }
+
+        public float Dq { get; set; }
+
+        public float Q { get; set; }
+
+        public float LastQ { get; set; }
+
+        public float Qmin { get; }
+
+        public float Qmax { get; }
+
+        public double Value { get; set; } // PSNR or size
+
+        public double LastValue { get; set; }
+
+        public double Target { get; }
+
+        public bool DoSizeSearch { get; }
+
+        public float ComputeNextQ()
+        {
+            float dq;
+            if (this.IsFirst)
+            {
+                dq = (this.Value > this.Target) ? -this.Dq : this.Dq;
+                this.IsFirst = false;
+            }
+            else if (this.Value != this.LastValue)
+            {
+                double slope = (this.Target - this.Value) / (this.LastValue - this.Value);
+                dq = (float)(slope * (this.LastQ - this.Q));
+            }
+            else
+            {
+                dq = 0.0f;  // we're done?!
+            }
+
+            // Limit variable to avoid large swings.
+            this.Dq = dq.Clamp(-30.0f, 30.0f);
+            this.LastQ = this.Q;
+            this.LastValue = this.Value;
+            this.Q = (this.Q + this.Dq).Clamp(this.Qmin, this.Qmax);
+
+            return this.Q;
+        }
+    }
+}

src/ImageSharp/Formats/WebP/Lossy/Vp8EncProba.cs

@@ -12,6 +12,11 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
         /// </summary>
         private const int MaxVariableLevel = 67;
 
+        /// <summary>
+        /// Value below which using skipProba is OK.
+        /// </summary>
+        private const int SkipProbaThreshold = 250;
+
         /// <summary>
         /// Initializes a new instance of the <see cref="Vp8EncProba"/> class.
         /// </summary>
@@ -30,6 +35,16 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
                 }
             }
 
+            this.Stats = new Vp8Stats[WebPConstants.NumTypes][];
+            for (int i = 0; i < this.Coeffs.Length; i++)
+            {
+                this.Stats[i] = new Vp8Stats[WebPConstants.NumBands];
+                for (int j = 0; j < this.Stats[i].Length; j++)
+                {
+                    this.Stats[i][j] = new Vp8Stats();
+                }
+            }
+
             this.LevelCost = new Vp8CostArray[WebPConstants.NumTypes][];
             for (int i = 0; i < this.LevelCost.Length; i++)
             {
@@ -74,17 +89,19 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
         public byte[] Segments { get; }
 
         /// <summary>
-        /// Gets the final probability of being skipped.
+        /// Gets or sets the final probability of being skipped.
         /// </summary>
-        public byte SkipProba { get; }
+        public byte SkipProba { get; set; }
 
         /// <summary>
-        /// Gets a value indicating whether to use the skip probability. Note: we always use SkipProba for now.
+        /// Gets or sets a value indicating whether to use the skip probability. Note: we always use SkipProba for now.
         /// </summary>
-        public bool UseSkipProba { get; }
+        public bool UseSkipProba { get; set; }
 
         public Vp8BandProbas[][] Coeffs { get; }
 
+        public Vp8Stats[][] Stats { get; }
+
         public Vp8CostArray[][] LevelCost { get; }
 
         public Vp8CostArray[][] RemappedCosts { get; }
@@ -132,7 +149,7 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
                     for (int ctx = 0; ctx < WebPConstants.NumCtx; ++ctx)
                     {
                         Span<ushort> dst = this.RemappedCosts[ctype][n].Costs.AsSpan(ctx * MaxVariableLevel, MaxVariableLevel);
-                        Span<ushort> src = this.LevelCost[ctype][WebPConstants.Bands[n]].Costs.AsSpan(ctx * MaxVariableLevel, MaxVariableLevel);
+                        Span<ushort> src = this.LevelCost[ctype][WebPConstants.Vp8EncBands[n]].Costs.AsSpan(ctx * MaxVariableLevel, MaxVariableLevel);
                         src.CopyTo(dst);
                     }
                 }
@@ -141,6 +158,87 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             this.Dirty = false;
         }
 
+        public int FinalizeTokenProbas()
+        {
+            bool hasChanged = false;
+            int size = 0;
+            for (int t = 0; t < WebPConstants.NumTypes; ++t)
+            {
+                for (int b = 0; b < WebPConstants.NumBands; ++b)
+                {
+                    for (int c = 0; c < WebPConstants.NumCtx; ++c)
+                    {
+                        for (int p = 0; p < WebPConstants.NumProbas; ++p)
+                        {
+                            var stats = this.Stats[t][b].Stats[c].Stats[p];
+                            int nb = (int)((stats >> 0) & 0xffff);
+                            int total = (int)((stats >> 16) & 0xffff);
+                            int updateProba = WebPLookupTables.CoeffsUpdateProba[t, b, c, p];
+                            int oldP = WebPLookupTables.DefaultCoeffsProba[t, b, c, p];
+                            int newP = this.CalcTokenProba(nb, total);
+                            int oldCost = this.BranchCost(nb, total, oldP) + this.BitCost(0, (byte)updateProba);
+                            int newCost = this.BranchCost(nb, total, newP) + this.BitCost(1, (byte)updateProba) + (8 * 256);
+                            bool useNewP = oldCost > newCost;
+                            size += this.BitCost(useNewP ? 1 : 0, (byte)updateProba);
+                            if (useNewP)
+                            {
+                                // Only use proba that seem meaningful enough.
+                                this.Coeffs[t][b].Probabilities[c].Probabilities[p] = (byte)newP;
+                                hasChanged |= newP != oldP;
+                                size += 8 * 256;
+                            }
+                            else
+                            {
+                                this.Coeffs[t][b].Probabilities[c].Probabilities[p] = (byte)oldP;
+                            }
+                        }
+                    }
+                }
+            }
+
+            this.Dirty = hasChanged;
+            return size;
+        }
+
+        public int FinalizeSkipProba(int mbw, int mbh)
+        {
+            int nbMbs = mbw * mbh;
+            int nbEvents = this.NbSkip;
+            this.SkipProba = (byte)this.CalcSkipProba(nbEvents, nbMbs);
+            this.UseSkipProba = this.SkipProba < SkipProbaThreshold;
+
+            int size = 256;
+            if (this.UseSkipProba)
+            {
+                size += (nbEvents * this.BitCost(1, this.SkipProba)) + ((nbMbs - nbEvents) * this.BitCost(0, this.SkipProba));
+                size += 8 * 256;   // cost of signaling the skipProba itself.
+            }
+
+            return size;
+        }
+
+        public void ResetTokenStats()
+        {
+            for (int t = 0; t < WebPConstants.NumTypes; ++t)
+            {
+                for (int b = 0; b < WebPConstants.NumBands; ++b)
+                {
+                    for (int c = 0; c < WebPConstants.NumCtx; ++c)
+                    {
+                        for (int p = 0; p < WebPConstants.NumProbas; ++p)
+                        {
+                            this.Stats[t][b].Stats[c].Stats[p] = 0;
+                        }
+                    }
+                }
+            }
+        }
+
+        private int CalcSkipProba(long nb, long total)
+        {
+            return (int)(total != 0 ? (total - nb) * 255 / total : 255);
+        }
+
         private int VariableLevelCost(int level, Span<byte> probas)
         {
             int pattern = WebPLookupTables.Vp8LevelCodes[level - 1][0];
@@ -160,6 +258,19 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             return cost;
         }
 
+        // Collect statistics and deduce probabilities for next coding pass.
+        // Return the total bit-cost for coding the probability updates.
+        private int CalcTokenProba(int nb, int total)
+        {
+            return nb != 0 ? (255 - (nb * 255 / total)) : 255;
+        }
+
+        // Cost of coding 'nb' 1's and 'total-nb' 0's using 'proba' probability.
+        private int BranchCost(int nb, int total, int proba)
+        {
+            return (nb * this.BitCost(1, (byte)proba)) + ((total - nb) * this.BitCost(0, (byte)proba));
+        }
+
         // Cost of coding one event with probability 'proba'.
         private int BitCost(int bit, byte proba)
         {

src/ImageSharp/Formats/WebP/Lossy/Vp8Encoder.cs

@@ -3,7 +3,6 @@
 
 using System;
 using System.Buffers;
-using System.Buffers.Binary;
 using System.IO;
 using System.Runtime.CompilerServices;
 
@@ -38,20 +37,25 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
         /// </summary>
         private readonly int method;
 
+        /// <summary>
+        /// Number of entropy-analysis passes (in [1..10]).
+        /// </summary>
+        private readonly int entropyPasses;
+
         /// <summary>
         /// Stride of the prediction plane (=4*mb_w + 1)
         /// </summary>
-        private int predsWidth;
+        private readonly int predsWidth;
 
         /// <summary>
         /// Macroblock width.
         /// </summary>
-        private int mbw;
+        private readonly int mbw;
 
         /// <summary>
         /// Macroblock height.
         /// </summary>
-        private int mbh;
+        private readonly int mbh;
 
         /// <summary>
         /// The segment features.
@@ -61,17 +65,21 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
         /// <summary>
         /// Contextual macroblock infos.
         /// </summary>
-        private Vp8MacroBlockInfo[] mbInfo;
+        private readonly Vp8MacroBlockInfo[] mbInfo;
 
         /// <summary>
         /// Probabilities.
         /// </summary>
-        private Vp8EncProba proba;
+        private readonly Vp8EncProba proba;
+
+        private readonly Vp8RdLevel rdOptLevel;
 
         private int dqUvDc;
 
         private int dqUvAc;
 
+        private int maxI4HeaderBits;
+
         /// <summary>
         /// Global susceptibility.
         /// </summary>
@@ -103,6 +111,17 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
 
         private const int MaxItersKMeans = 6;
 
+        // Convergence is considered reached if dq < DqLimit
+        private const float DqLimit = 0.4f;
+
+        private const ulong Partition0SizeLimit = (WebPConstants.Vp8MaxPartition0Size - 2048UL) << 11;
+
+        private const long HeaderSizeEstimate = WebPConstants.RiffHeaderSize + WebPConstants.ChunkHeaderSize + WebPConstants.Vp8FrameHeaderSize;
+
+        private const int QMin = 0;
+
+        private const int QMax = 100;
+
         /// <summary>
         /// Initializes a new instance of the <see cref="Vp8Encoder"/> class.
         /// </summary>
@@ -111,11 +130,17 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
         /// <param name="height">The height of the input image.</param>
         /// <param name="quality">The encoding quality.</param>
         /// <param name="method">Quality/speed trade-off (0=fast, 6=slower-better).</param>
-        public Vp8Encoder(MemoryAllocator memoryAllocator, int width, int height, int quality, int method)
+        /// <param name="entropyPasses">Number of entropy-analysis passes (in [1..10]).</param>
+        public Vp8Encoder(MemoryAllocator memoryAllocator, int width, int height, int quality, int method, int entropyPasses)
         {
             this.memoryAllocator = memoryAllocator;
             this.quality = quality.Clamp(0, 100);
             this.method = method.Clamp(0, 6);
+            this.entropyPasses = entropyPasses.Clamp(1, 10);
+            this.rdOptLevel = (method >= 6) ? Vp8RdLevel.RdOptTrellisAll
+                : (method >= 5) ? Vp8RdLevel.RdOptTrellis
+                : (method >= 3) ? Vp8RdLevel.RdOptBasic
+                : Vp8RdLevel.RdOptNone;
 
             var pixelCount = width * height;
             this.mbw = (width + 15) >> 4;
@@ -130,6 +155,12 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             this.MbHeaderLimit = 256 * 510 * 8 * 1024 / (this.mbw * this.mbh);
             int predSize = (((4 * this.mbw) + 1) * ((4 * this.mbh) + 1)) + this.predsWidth + 1;
 
+            // TODO: make partition_limit configurable?
+            int limit = 100; // original code: limit = 100 - config->partition_limit;
+            this.maxI4HeaderBits =
+                256 * 16 * 16 * // upper bound: up to 16bit per 4x4 block
+                (limit * limit) / (100 * 100);  // ... modulated with a quadratic curve.
+
             this.mbInfo = new Vp8MacroBlockInfo[this.mbw * this.mbh];
             for (int i = 0; i < this.mbInfo.Length; i++)
             {
@@ -225,20 +256,22 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             this.alpha = this.MacroBlockAnalysis(width, height, it, y, u, v, yStride, uvStride, alphas, out this.uvAlpha);
 
             // Analysis is done, proceed to actual encoding.
-
-            // TODO: EncodeAlpha();
             this.segmentHeader = new Vp8EncSegmentHeader(4);
             this.AssignSegments(segmentInfos, alphas);
-            this.SetSegmentParams(segmentInfos);
+            this.SetSegmentParams(segmentInfos, this.quality);
             this.SetSegmentProbas(segmentInfos);
             this.ResetStats();
+
+            // TODO: EncodeAlpha();
+            // Stats-collection loop.
+            this.StatLoop(width, height, yStride, uvStride, segmentInfos);
             it.Init();
             it.InitFilter();
             do
             {
                 var info = new Vp8ModeScore();
                 it.Import(y, u, v, yStride, uvStride, width, height);
-                if (!this.Decimate(it, segmentInfos, info))
+                if (!this.Decimate(it, segmentInfos, info, this.rdOptLevel))
                 {
                     this.CodeResiduals(it, info);
                 }
@@ -266,6 +299,139 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             this.Preds.Dispose();
         }
 
+        /// <summary>
+        /// Only collect statistics(number of skips, token usage, ...).
+        /// This is used for deciding optimal probabilities. It also modifies the
+        /// quantizer value if some target (size, PSNR) was specified.
+        /// </summary>
+        private void StatLoop(int width, int height, int yStride, int uvStride, Vp8SegmentInfo[] segmentInfos)
+        {
+            int targetSize = 0; // TODO: target size is hardcoded.
+            float targetPsnr = 0.0f; // TDOO: targetPsnr is hardcoded.
+            int method = this.method;
+            bool doSearch = false; // TODO: doSearch hardcoded for now.
+            bool fastProbe = (method == 0 || method == 3) && !doSearch;
+            int numPassLeft = this.entropyPasses;
+            Vp8RdLevel rdOpt = (method >= 3 || doSearch) ? Vp8RdLevel.RdOptBasic : Vp8RdLevel.RdOptNone;
+            int nbMbs = this.mbw * this.mbh;
+
+            var stats = new PassStats(targetSize, targetPsnr, QMin, QMax, this.quality);
+            this.proba.ResetTokenStats();
+
+            // Fast mode: quick analysis pass over few mbs. Better than nothing.
+            if (fastProbe)
+            {
+                if (method == 3)
+                {
+                    // We need more stats for method 3 to be reliable.
+                    nbMbs = (nbMbs > 200) ? nbMbs >> 1 : 100;
+                }
+                else
+                {
+                    nbMbs = (nbMbs > 200) ? nbMbs >> 2 : 50;
+                }
+            }
+
+            while (numPassLeft-- > 0)
+            {
+                bool isLastPass = (MathF.Abs(stats.Dq) <= DqLimit) || (numPassLeft == 0) || (this.maxI4HeaderBits == 0);
+                var sizeP0 = this.OneStatPass(width, height, yStride, uvStride, rdOpt, nbMbs, stats, segmentInfos);
+                if (sizeP0 == 0)
+                {
+                    return;
+                }
+
+                if (this.maxI4HeaderBits > 0 && sizeP0 > (long)Partition0SizeLimit)
+                {
+                    ++numPassLeft;
+                    this.maxI4HeaderBits >>= 1;  // strengthen header bit limitation...
+                    continue;                        // ...and start over
+                }
+
+                if (isLastPass)
+                {
+                    break;
+                }
+
+                // If no target size: just do several pass without changing 'q'
+                if (doSearch)
+                {
+                    stats.ComputeNextQ();
+                    if (MathF.Abs(stats.Dq) <= DqLimit)
+                    {
+                        break;
+                    }
+                }
+            }
+
+            if (!doSearch || !stats.DoSizeSearch)
+            {
+                // Need to finalize probas now, since it wasn't done during the search.
+                this.proba.FinalizeSkipProba(this.mbw, this.mbh);
+                this.proba.FinalizeTokenProbas();
+            }
+
+            this.proba.CalculateLevelCosts();  // finalize costs
+        }
+
+        private long OneStatPass(int width, int height, int yStride, int uvStride, Vp8RdLevel rdOpt, int nbMbs, PassStats stats, Vp8SegmentInfo[] segmentInfos)
+        {
+            Span<byte> y = this.Y.GetSpan();
+            Span<byte> u = this.U.GetSpan();
+            Span<byte> v = this.V.GetSpan();
+            var it = new Vp8EncIterator(this.YTop, this.UvTop, this.Preds, this.Nz, this.mbInfo, this.mbw, this.mbh);
+            long size = 0;
+            long sizeP0 = 0;
+            long distortion = 0;
+            long pixelCount = nbMbs * 384;
+
+            this.SetLoopParams(segmentInfos, stats.Q);
+            do
+            {
+                var info = new Vp8ModeScore();
+                it.Import(y, u, v, yStride, uvStride, width, height);
+                if (this.Decimate(it, segmentInfos, info, rdOpt))
+                {
+                    // Just record the number of skips and act like skipProba is not used.
+                    ++this.proba.NbSkip;
+                }
+
+                this.RecordResiduals(it, info);
+                size += info.R + info.H;
+                sizeP0 += info.H;
+                distortion += info.D;
+
+                it.SaveBoundary();
+            }
+            while (it.Next());
+
+            sizeP0 += this.segmentHeader.Size;
+            if (stats.DoSizeSearch)
+            {
+                size += this.proba.FinalizeSkipProba(this.mbw, this.mbh);
+                size += this.proba.FinalizeTokenProbas();
+                size = ((size + sizeP0 + 1024) >> 11) + HeaderSizeEstimate;
+                stats.Value = size;
+            }
+            else
+            {
+                stats.Value = this.GetPsnr(distortion, pixelCount);
+            }
+
+            return sizeP0;
+        }
+
+        private void SetLoopParams(Vp8SegmentInfo[] dqm, float q)
+        {
+            // Setup segment quantizations and filters.
+            this.SetSegmentParams(dqm, q);
+
+            // Compute segment probabilities.
+            this.SetSegmentProbas(dqm);
+
+            this.ResetStats();
+        }
+
         private void ResetBoundaryPredictions()
         {
             Span<byte> top = this.Preds.GetSpan();
@@ -416,12 +582,12 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             }
         }
 
-        private void SetSegmentParams(Vp8SegmentInfo[] dqm)
+        private void SetSegmentParams(Vp8SegmentInfo[] dqm, float quality)
         {
             int nb = this.segmentHeader.NumSegments;
             int snsStrength = 50; // TODO: Spatial Noise Shaping, hardcoded for now.
             double amp = WebPConstants.SnsToDq * snsStrength / 100.0d / 128.0d;
-            double cBase = this.QualityToCompression(this.quality / 100.0d);
+            double cBase = this.QualityToCompression(quality / 100.0d);
             for (int i = 0; i < nb; ++i)
             {
                 // We modulate the base coefficient to accommodate for the quantization
@@ -488,6 +654,8 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
                 m.Uv.Q[1] = WebPLookupTables.AcTable[this.Clip(q + this.dqUvAc, 0, 127)];
 
                 var qi4 = m.Y1.Expand(0);
+                var qi16 = m.Y2.Expand(1);
+                var quv = m.Uv.Expand(2);
 
                 m.I4Penalty = 1000 * qi4 * qi4;
             }
@@ -541,7 +709,7 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             return bestAlpha; // Mixed susceptibility (not just luma).
         }
 
-        private bool Decimate(Vp8EncIterator it, Vp8SegmentInfo[] segmentInfos, Vp8ModeScore rd)
+        private bool Decimate(Vp8EncIterator it, Vp8SegmentInfo[] segmentInfos, Vp8ModeScore rd, Vp8RdLevel rdOpt)
         {
             rd.InitScore();
 
@@ -730,7 +898,7 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
                 for (x = 0; x < 4; ++x)
                 {
                     int ctx = it.TopNz[x] + it.LeftNz[y];
-                    residual.SetCoeffs(rd.YAcLevels.AsSpan(x + (y * 4), 16));
+                    residual.SetCoeffs(rd.YAcLevels.AsSpan(16 * (x + (y * 4)), 16));
                     int res = this.bitWriter.PutCoeffs(ctx, residual);
                     it.TopNz[x] = it.LeftNz[y] = res;
                 }
@@ -747,7 +915,7 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
                     for (x = 0; x < 2; ++x)
                     {
                         int ctx = it.TopNz[4 + ch + x] + it.LeftNz[4 + ch + y];
-                        residual.SetCoeffs(rd.UvLevels.AsSpan((ch * 2) + x + (y * 2), 16));
+                        residual.SetCoeffs(rd.UvLevels.AsSpan(16 * ((ch * 2) + x + (y * 2)), 16));
                         var res = this.bitWriter.PutCoeffs(ctx, residual);
                         it.TopNz[4 + ch + x] = it.LeftNz[4 + ch + y] = res;
                     }
@@ -762,6 +930,67 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             it.BytesToNz();
         }
 
+        /// <summary>
+        /// Same as CodeResiduals, but doesn't actually write anything.
+        /// Instead, it just records the event distribution.
+        /// </summary>
+        private void RecordResiduals(Vp8EncIterator it, Vp8ModeScore rd)
+        {
+            int x, y, ch;
+            var residual = new Vp8Residual();
+            bool i16 = it.CurrentMacroBlockInfo.MacroBlockType == Vp8MacroBlockType.I16X16;
+            int segment = it.CurrentMacroBlockInfo.Segment;
+
+            it.NzToBytes();
+
+            if (i16)
+            {
+                // i16x16
+                residual.Init(0, 1, this.proba);
+                residual.SetCoeffs(rd.YDcLevels);
+                var res = residual.RecordCoeffs(it.TopNz[8] + it.LeftNz[8]);
+                it.TopNz[8] = res;
+                it.LeftNz[8] = res;
+                residual.Init(1, 0, this.proba);
+            }
+            else
+            {
+                residual.Init(0, 3, this.proba);
+            }
+
+            // luma-AC
+            for (y = 0; y < 4; ++y)
+            {
+                for (x = 0; x < 4; ++x)
+                {
+                    int ctx = it.TopNz[x] + it.LeftNz[y];
+                    residual.SetCoeffs(rd.YAcLevels.AsSpan(16 * (x + (y * 4)), 16));
+                    var res = residual.RecordCoeffs(ctx);
+                    it.TopNz[x] = res;
+                    it.LeftNz[y] = res;
+                }
+            }
+
+            // U/V
+            residual.Init(0, 2, this.proba);
+            for (ch = 0; ch <= 2; ch += 2)
+            {
+                for (y = 0; y < 2; ++y)
+                {
+                    for (x = 0; x < 2; ++x)
+                    {
+                        int ctx = it.TopNz[4 + ch + x] + it.LeftNz[4 + ch + y];
+                        residual.SetCoeffs(rd.UvLevels.AsSpan(16 * ((ch * 2) + x + (y * 2)), 16));
+                        var res = residual.RecordCoeffs(ctx);
+                        it.TopNz[4 + ch + x] = res;
+                        it.LeftNz[4 + ch + y] = res;
+                    }
+                }
+            }
+
+            it.BytesToNz();
+        }
+
         private int ReconstructIntra16(Vp8EncIterator it, Vp8SegmentInfo dqm, Vp8ModeScore rd, Span<byte> yuvOut, int mode)
         {
             Span<byte> reference = it.YuvP.AsSpan(Vp8EncIterator.Vp8I16ModeOffsets[mode]);
@@ -785,7 +1014,7 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
                 // Zero-out the first coeff, so that: a) nz is correct below, and
                 // b) finding 'last' non-zero coeffs in SetResidualCoeffs() is simplified.
                 tmp[n * 16] = tmp[(n + 1) * 16] = 0;
-                nz |= this.Quantize2Blocks(tmpSpan.Slice(n * 16), rd.YAcLevels.AsSpan(n, 32), dqm.Y1) << n;
+                nz |= this.Quantize2Blocks(tmpSpan.Slice(n * 16), rd.YAcLevels.AsSpan(n * 16, 32), dqm.Y1) << n;
             }
 
             // Transform back.
@@ -1400,6 +1629,12 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
             return v;
         }
 
+        [MethodImpl(InliningOptions.ShortMethod)]
+        private double GetPsnr(long mse, long size)
+        {
+            return (mse > 0 && size > 0) ? 10.0f * Math.Log10(255.0f * 255.0f * size / mse) : 99;
+        }
+
         [MethodImpl(InliningOptions.ShortMethod)]
         private int QuantDiv(uint n, uint iQ, uint b)
         {

src/ImageSharp/Formats/WebP/Lossy/Vp8ModeScore.cs

@@ -17,7 +17,7 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
         {
             this.YDcLevels = new short[16];
             this.YAcLevels = new short[16 * 16];
-            this.UvLevels = new short[4 + (4 * 16)];
+            this.UvLevels = new short[(4 + 4) * 16];
 
             this.ModesI4 = new byte[16];
         }

src/ImageSharp/Formats/WebP/Lossy/Vp8RDLevel.cs

@@ -0,0 +1,31 @@
+// Copyright (c) Six Labors.
+// Licensed under the Apache License, Version 2.0.
+
+namespace SixLabors.ImageSharp.Formats.WebP.Lossy
+{
+    /// <summary>
+    /// Rate-distortion optimization levels
+    /// </summary>
+    internal enum Vp8RdLevel
+    {
+        /// <summary>
+        /// No rd-opt.
+        /// </summary>
+        RdOptNone = 0,
+
+        /// <summary>
+        /// Basic scoring (no trellis).
+        /// </summary>
+        RdOptBasic = 1,
+
+        /// <summary>
+        /// Perform trellis-quant on the final decision only.
+        /// </summary>
+        RdOptTrellis = 2,
+
+        /// <summary>
+        /// Trellis-quant for every scoring (much slower).
+        /// </summary>
+        RdOptTrellisAll = 3
+    }
+}

src/ImageSharp/Formats/WebP/Lossy/Vp8Residual.cs

@@ -10,6 +10,8 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
     /// </summary>
     internal class Vp8Residual
     {
+        private const int MaxVariableLevel = 67;
+
         public int First { get; set; }
 
         public int Last { get; set; }
@@ -20,14 +22,16 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
 
         public Vp8BandProbas[] Prob { get; set; }
 
+        public Vp8Stats[] Stats { get; set; }
+
         public void Init(int first, int coeffType, Vp8EncProba prob)
         {
             this.First = first;
             this.CoeffType = coeffType;
             this.Prob = prob.Coeffs[this.CoeffType];
+            this.Stats = prob.Stats[this.CoeffType];
 
             // TODO:
-            // res->stats = enc->proba_.stats_[coeff_type];
             // res->costs = enc->proba_.remapped_costs_[coeff_type];
         }
 
@@ -46,5 +50,79 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
 
             this.Coeffs = coeffs.ToArray();
         }
+
+        // Simulate block coding, but only record statistics.
+        // Note: no need to record the fixed probas.
+        public int RecordCoeffs(int ctx)
+        {
+            int n = this.First;
+            Vp8StatsArray s = this.Stats[n].Stats[ctx];
+            if (this.Last < 0)
+            {
+                this.RecordStats(0, s,  0);
+                return 0;
+            }
+
+            while (n <= this.Last)
+            {
+                int v;
+                this.RecordStats(1, s, 0);  // order of record doesn't matter
+                while ((v = this.Coeffs[n++]) == 0)
+                {
+                    this.RecordStats(0, s, 1);
+                    s = this.Stats[WebPConstants.Vp8EncBands[n]].Stats[0];
+                    this.RecordStats(1, s, 1);
+                    if (this.RecordStats((v + 1) > 2u ? 1 : 0, s, 2) == 0)
+                    {
+                        // v = -1 or 1
+                        s = this.Stats[WebPConstants.Vp8EncBands[n]].Stats[1];
+                    }
+                    else
+                    {
+                        v = Math.Abs(v);
+                        if (v > MaxVariableLevel)
+                        {
+                            v = MaxVariableLevel;
+                        }
+
+                        int bits = WebPLookupTables.Vp8LevelCodes[v - 1][1];
+                        int pattern = WebPLookupTables.Vp8LevelCodes[v - 1][0];
+                        int i;
+                        for (i = 0; (pattern >>= 1) != 0; ++i)
+                        {
+                            int mask = 2 << i;
+                            if ((pattern & 1) != 0)
+                            {
+                                this.RecordStats(bits & mask, s, 3 + i);
+                            }
+                        }
+
+                        s = this.Stats[WebPConstants.Vp8EncBands[n]].Stats[2];
+                    }
+                }
+            }
+
+            if (n < 16)
+            {
+                this.RecordStats(0, s,  0);
+            }
+
+            return 1;
+        }
+
+        private int RecordStats(int bit, Vp8StatsArray statsArr, int idx)
+        {
+            // An overflow is inbound. Note we handle this at 0xfffe0000u instead of
+            // 0xffff0000u to make sure p + 1u does not overflow.
+            if (statsArr.Stats[idx] >= 0xfffe0000u)
+            {
+                statsArr.Stats[idx] = ((statsArr.Stats[idx] + 1u) >> 1) & 0x7fff7fffu;  // -> divide the stats by 2.
+            }
+
+            // record bit count (lower 16 bits) and increment total count (upper 16 bits).
+            statsArr.Stats[idx] += 0x00010000u + (uint)bit;
+
+            return bit;
+        }
     }
 }

src/ImageSharp/Formats/WebP/Lossy/Vp8Stats.cs

@@ -0,0 +1,22 @@
+// Copyright (c) Six Labors.
+// Licensed under the Apache License, Version 2.0.
+
+namespace SixLabors.ImageSharp.Formats.WebP.Lossy
+{
+    internal class Vp8Stats
+    {
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Vp8Stats"/> class.
+        /// </summary>
+        public Vp8Stats()
+        {
+            this.Stats = new Vp8StatsArray[WebPConstants.NumCtx];
+            for (int i = 0; i < WebPConstants.NumCtx; i++)
+            {
+                this.Stats[i] = new Vp8StatsArray();
+            }
+        }
+
+        public Vp8StatsArray[] Stats { get; }
+    }
+}

src/ImageSharp/Formats/WebP/Lossy/Vp8StatsArray.cs

@@ -0,0 +1,18 @@
+// Copyright (c) Six Labors.
+// Licensed under the Apache License, Version 2.0.
+
+namespace SixLabors.ImageSharp.Formats.WebP.Lossy
+{
+    internal class Vp8StatsArray
+    {
+        /// <summary>
+        /// Initializes a new instance of the <see cref="Vp8StatsArray"/> class.
+        /// </summary>
+        public Vp8StatsArray()
+        {
+            this.Stats = new uint[WebPConstants.NumProbas];
+        }
+
+        public uint[] Stats { get; }
+    }
+}

src/ImageSharp/Formats/WebP/Lossy/WebPLossyDecoder.cs

@@ -1278,7 +1278,7 @@ namespace SixLabors.ImageSharp.Formats.WebP.Lossy
 
                 for (int b = 0; b < 16 + 1; ++b)
                 {
-                    proba.BandsPtr[t][b] = proba.Bands[t, WebPConstants.Bands[b]];
+                    proba.BandsPtr[t][b] = proba.Bands[t, WebPConstants.Vp8EncBands[b]];
                 }
             }
 

src/ImageSharp/Formats/WebP/WebPConstants.cs

@@ -6,7 +6,7 @@ using System.Collections.Generic;
 namespace SixLabors.ImageSharp.Formats.WebP
 {
     /// <summary>
-    /// Constants used for decoding VP8 and VP8L bitstreams.
+    /// Constants used for encoding and decoding VP8 and VP8L bitstreams.
     /// </summary>
     internal static class WebPConstants
     {
@@ -78,11 +78,21 @@ namespace SixLabors.ImageSharp.Formats.WebP
         /// </summary>
         public const int Vp8LImageSizeBits = 14;
 
+        /// <summary>
+        /// Size of the frame header within VP8 data.
+        /// </summary>
+        public const int Vp8FrameHeaderSize = 10;
+
         /// <summary>
         /// Size of a chunk header.
         /// </summary>
         public const int ChunkHeaderSize = 8;
 
+        /// <summary>
+        /// Size of the RIFF header ("RIFFnnnnWEBP").
+        /// </summary>
+        public const int RiffHeaderSize = 12;
+
         /// <summary>
         /// Size of a chunk tag (e.g. "VP8L").
         /// </summary>
@@ -241,6 +251,11 @@ namespace SixLabors.ImageSharp.Formats.WebP
 
         public const int QFix = 17;
 
+        /// <summary>
+        /// Max size of mode partition.
+        /// </summary>
+        public const int Vp8MaxPartition0Size = 1 << 19;
+
         public static readonly short[] Vp8FixedCostsUv = { 302, 984, 439, 642 };
 
         public static readonly short[] Vp8FixedCostsI16 = { 663, 919, 872, 919 };
@@ -258,7 +273,7 @@ namespace SixLabors.ImageSharp.Formats.WebP
         public static readonly byte[] FilterExtraRows = { 0, 2, 8 };
 
         // Paragraph 9.9
-        public static readonly int[] Bands =
+        public static readonly int[] Vp8EncBands =
         {
             0, 1, 2, 3, 6, 4, 5, 6, 6, 6, 6, 6, 6, 6, 6, 7, 0
         };

src/ImageSharp/Formats/WebP/WebPEncoderCore.cs

@@ -48,6 +48,11 @@ namespace SixLabors.ImageSharp.Formats.WebP
         /// </summary>
         private readonly int method;
 
+        /// <summary>
+        /// The number of entropy-analysis passes (in [1..10]).
+        /// </summary>
+        private readonly int entropyPasses;
+
         /// <summary>
         /// Initializes a new instance of the <see cref="WebPEncoderCore"/> class.
         /// </summary>
@@ -60,6 +65,7 @@ namespace SixLabors.ImageSharp.Formats.WebP
             this.lossy = options.Lossy;
             this.quality = options.Quality;
             this.method = options.Method;
+            this.entropyPasses = options.EntropyPasses;
         }
 
         /// <summary>
@@ -79,7 +85,7 @@ namespace SixLabors.ImageSharp.Formats.WebP
 
             if (this.lossy)
             {
-                var enc = new Vp8Encoder(this.memoryAllocator, image.Width, image.Height, this.quality, this.method);
+                var enc = new Vp8Encoder(this.memoryAllocator, image.Width, image.Height, this.quality, this.method, this.entropyPasses);
                 enc.Encode(image, stream);
             }
             else