< Summary

Information
Class: JeremyAnsel.ColorQuant.WuColorQuantizer
Assembly: JeremyAnsel.ColorQuant
File(s): File 1: https://raw.githubusercontent.com/JeremyAnsel/JeremyAnsel.ColorQuant/6d79217e72af9e3af1a8a29c606732adac1e8d87/JeremyAnsel.ColorQuant/JeremyAnsel.ColorQuant/WuColorQuantizer.cs
File 2: https://raw.githubusercontent.com/JeremyAnsel/JeremyAnsel.ColorQuant/6d79217e72af9e3af1a8a29c606732adac1e8d87/JeremyAnsel.ColorQuant/JeremyAnsel.ColorQuant/WuColorQuantizer2.cs
Line coverage
96%
Covered lines: 397
Uncovered lines: 14
Coverable lines: 411
Total lines: 1162
Line coverage: 96.5%
Branch coverage
93%
Covered branches: 134
Total branches: 144
Branch coverage: 93%
Method coverage

Feature is only available for sponsors

Upgrade to PRO version

Metrics

MethodBranch coverage Crap Score Cyclomatic complexity Line coverage
File 1: .ctor()100%11100%
File 1: Quantize(...)100%11100%
File 1: Quantize(...)100%66100%
File 1: GetIndex(...)100%11100%
File 1: Volume(...)100%11100%
File 1: Bottom(...)75%4492.3%
File 1: Top(...)75%4492.3%
File 1: Clear()100%11100%
File 1: Build3DHistogram(...)100%22100%
File 1: Get3DMoments()100%66100%
File 1: Variance(...)100%11100%
File 1: Maximize(...)100%88100%
File 1: Cut(...)100%1414100%
File 1: Mark(...)100%66100%
File 1: BuildCube(...)87.5%161695.83%
File 1: GenerateResult(...)83.33%6680%
File 2: .ctor()100%11100%
File 2: Quantize(...)100%11100%
File 2: Quantize(...)100%66100%
File 2: GetIndex(...)100%11100%
File 2: Volume(...)100%11100%
File 2: Bottom(...)75%4492.3%
File 2: Top(...)75%4492.3%
File 2: Clear()100%11100%
File 2: Build3DHistogram(...)100%22100%
File 2: Get3DMoments()100%66100%
File 2: Variance(...)100%11100%
File 2: Maximize(...)100%88100%
File 2: Cut(...)100%1414100%
File 2: Mark(...)100%66100%
File 2: BuildCube(...)87.5%161695.83%
File 2: GenerateResult(...)83.33%6680.95%

File(s)

https://raw.githubusercontent.com/JeremyAnsel/JeremyAnsel.ColorQuant/6d79217e72af9e3af1a8a29c606732adac1e8d87/JeremyAnsel.ColorQuant/JeremyAnsel.ColorQuant/WuColorQuantizer.cs

#LineLine coverage
 1// <copyright file="WuColorQuantizer.cs" company="Jérémy Ansel">
 2// Copyright (c) 2014-2019 Jérémy Ansel
 3// </copyright>
 4// <license>
 5// Licensed under the MIT license. See LICENSE.txt
 6// </license>
 7
 8#if !NET8_0_OR_GREATER
 9namespace JeremyAnsel.ColorQuant
 10{
 11    using System;
 12    using System.Diagnostics.CodeAnalysis;
 13
 14    /// <summary>
 15    /// A Wu's color quantizer.
 16    /// </summary>
 17    /// <remarks>
 18    /// <para>
 19    /// Based on C Implementation of Xiaolin Wu's Color Quantizer (v. 2)
 20    /// (see Graphics Gems volume II, pages 126-133)
 21    /// (<see href="http://www.ece.mcmaster.ca/~xwu/cq.c"/>).
 22    /// </para>
 23    /// <para>
 24    /// Algorithm: Greedy orthogonal bipartition of RGB space for variance
 25    /// minimization aided by inclusion-exclusion tricks.
 26    /// For speed no nearest neighbor search is done. Slightly
 27    /// better performance can be expected by more sophisticated
 28    /// but more expensive versions.
 29    /// </para>
 30    /// </remarks>
 31    [SuppressMessage("Microsoft.Naming", "CA1709:IdentifiersShouldBeCasedCorrectly", MessageId = "Wu", Justification = "
 32    public sealed class WuColorQuantizer : IColorQuantizer
 33    {
 34        /// <summary>
 35        /// The maximum color count for the quantize.
 36        /// </summary>
 37        public const int MaxColors = 256;
 38
 39        /// <summary>
 40        /// The index bits.
 41        /// </summary>
 42        private const int IndexBits = 7;
 43
 44        /// <summary>
 45        /// The index count.
 46        /// </summary>
 47        private const int IndexCount = (1 << IndexBits) + 1;
 48
 49        /// <summary>
 50        /// The table length.
 51        /// </summary>
 52        private const int TableLength = IndexCount * IndexCount * IndexCount;
 53
 54        /// <summary>
 55        /// Moment of <c>P(c)</c>.
 56        /// </summary>
 157        private readonly long[] vwt = new long[TableLength];
 58
 59        /// <summary>
 60        /// Moment of <c>r*P(c)</c>.
 61        /// </summary>
 162        private readonly long[] vmr = new long[TableLength];
 63
 64        /// <summary>
 65        /// Moment of <c>g*P(c)</c>.
 66        /// </summary>
 167        private readonly long[] vmg = new long[TableLength];
 68
 69        /// <summary>
 70        /// Moment of <c>b*P(c)</c>.
 71        /// </summary>
 172        private readonly long[] vmb = new long[TableLength];
 73
 74        /// <summary>
 75        /// Moment of <c>c^2*P(c)</c>.
 76        /// </summary>
 177        private readonly double[] m2 = new double[TableLength];
 78
 79        /// <summary>
 80        /// Color space tag.
 81        /// </summary>
 182        private readonly byte[] tag = new byte[TableLength];
 83
 84        /// <summary>
 85        /// Quantizes an image.
 86        /// </summary>
 87        /// <param name="image">The image (XRGB).</param>
 88        /// <returns>The result.</returns>
 89        public ColorQuantizerResult Quantize(byte[] image)
 90        {
 1191            return this.Quantize(image, 256);
 92        }
 93
 94        /// <summary>
 95        /// Quantizes an image.
 96        /// </summary>
 97        /// <param name="image">The image (XRGB).</param>
 98        /// <param name="colorCount">The color count.</param>
 99        /// <returns>The result.</returns>
 100        public ColorQuantizerResult Quantize(byte[] image, int colorCount)
 101        {
 14102            if (image == null)
 103            {
 2104                throw new ArgumentNullException(nameof(image));
 105            }
 106
 12107            if (colorCount < 1 || colorCount > 256)
 108            {
 2109                throw new ArgumentOutOfRangeException(nameof(colorCount));
 110            }
 111
 10112            this.Clear();
 113
 10114            this.Build3DHistogram(image);
 10115            this.Get3DMoments();
 116
 117            Box[] cube;
 10118            this.BuildCube(out cube, ref colorCount);
 119
 10120            return this.GenerateResult(image, colorCount, cube);
 121        }
 122
 123        /// <summary>
 124        /// Gets an index.
 125        /// </summary>
 126        /// <param name="r">The red value.</param>
 127        /// <param name="g">The green value.</param>
 128        /// <param name="b">The blue value.</param>
 129        /// <returns>The index.</returns>
 130        private static int GetIndex(int r, int g, int b)
 131        {
 199588134132            return (r << (IndexBits * 2)) + (r << (IndexBits + 1)) + (g << IndexBits) + r + g + b;
 133        }
 134
 135        /// <summary>
 136        /// Computes sum over a box of any given statistic.
 137        /// </summary>
 138        /// <param name="cube">The cube.</param>
 139        /// <param name="moment">The moment.</param>
 140        /// <returns>The result.</returns>
 141        private static double Volume(Box cube, long[] moment)
 142        {
 24472143            return moment[cube.Index111]
 24472144               - moment[cube.Index110]
 24472145               - moment[cube.Index101]
 24472146               + moment[cube.Index100]
 24472147               - moment[cube.Index011]
 24472148               + moment[cube.Index010]
 24472149               + moment[cube.Index001]
 24472150               - moment[cube.Index000];
 151        }
 152
 153        /// <summary>
 154        /// Computes part of Volume(cube, moment) that doesn't depend on r1, g1, or b1 (depending on direction).
 155        /// </summary>
 156        /// <param name="cube">The cube.</param>
 157        /// <param name="direction">The direction.</param>
 158        /// <param name="moment">The moment.</param>
 159        /// <returns>The result.</returns>
 160        private static long Bottom(Box cube, int direction, long[] moment)
 161        {
 162            switch (direction)
 163            {
 164                // Red
 165                case 2:
 9184166                    return -moment[cube.Index011]
 9184167                        + moment[cube.Index010]
 9184168                        + moment[cube.Index001]
 9184169                        - moment[cube.Index000];
 170
 171                // Green
 172                case 1:
 9184173                    return -moment[cube.Index101]
 9184174                        + moment[cube.Index100]
 9184175                        + moment[cube.Index001]
 9184176                        - moment[cube.Index000];
 177
 178                // Blue
 179                case 0:
 9184180                    return -moment[cube.Index110]
 9184181                        + moment[cube.Index100]
 9184182                        + moment[cube.Index010]
 9184183                        - moment[cube.Index000];
 184
 185                default:
 0186                    throw new ArgumentOutOfRangeException(nameof(direction));
 187            }
 188        }
 189
 190        /// <summary>
 191        /// Computes remainder of Volume(cube, moment), substituting position for r1, g1, or b1 (depending on direction)
 192        /// </summary>
 193        /// <param name="cube">The cube.</param>
 194        /// <param name="direction">The direction.</param>
 195        /// <param name="position">The position.</param>
 196        /// <param name="moment">The moment.</param>
 197        /// <returns>The result.</returns>
 198        private static long Top(Box cube, int direction, int position, long[] moment)
 199        {
 200            switch (direction)
 201            {
 202                // Red
 203                case 2:
 159866204                    return moment[GetIndex(position, cube.G1, cube.B1)]
 159866205                       - moment[GetIndex(position, cube.G1, cube.B0)]
 159866206                       - moment[GetIndex(position, cube.G0, cube.B1)]
 159866207                       + moment[GetIndex(position, cube.G0, cube.B0)];
 208
 209                // Green
 210                case 1:
 215797211                    return moment[GetIndex(cube.R1, position, cube.B1)]
 215797212                       - moment[GetIndex(cube.R1, position, cube.B0)]
 215797213                       - moment[GetIndex(cube.R0, position, cube.B1)]
 215797214                       + moment[GetIndex(cube.R0, position, cube.B0)];
 215
 216                // Blue
 217                case 0:
 237658218                    return moment[GetIndex(cube.R1, cube.G1, position)]
 237658219                       - moment[GetIndex(cube.R1, cube.G0, position)]
 237658220                       - moment[GetIndex(cube.R0, cube.G1, position)]
 237658221                       + moment[GetIndex(cube.R0, cube.G0, position)];
 222
 223                default:
 0224                    throw new ArgumentOutOfRangeException(nameof(direction));
 225            }
 226        }
 227
 228        /// <summary>
 229        /// Clears the tables.
 230        /// </summary>
 231        private void Clear()
 232        {
 10233            this.vwt.AsSpan().Clear();
 10234            this.vmr.AsSpan().Clear();
 10235            this.vmg.AsSpan().Clear();
 10236            this.vmb.AsSpan().Clear();
 10237            this.m2.AsSpan().Clear();
 10238            this.tag.AsSpan().Clear();
 10239        }
 240
 241        /// <summary>
 242        /// Builds a 3-D color histogram of <c>counts, r/g/b, c^2</c>.
 243        /// </summary>
 244        /// <param name="image">The image.</param>
 245        private void Build3DHistogram(byte[] image)
 246        {
 134220310247            for (int i = 0; i < image.Length; i += 4)
 248            {
 67110145249                int r = image[i + 2];
 67110145250                int g = image[i + 1];
 67110145251                int b = image[i];
 252
 67110145253                int inr = r >> (8 - IndexBits);
 67110145254                int ing = g >> (8 - IndexBits);
 67110145255                int inb = b >> (8 - IndexBits);
 256
 67110145257                int ind = GetIndex(inr + 1, ing + 1, inb + 1);
 258
 67110145259                this.vwt[ind]++;
 67110145260                this.vmr[ind] += r;
 67110145261                this.vmg[ind] += g;
 67110145262                this.vmb[ind] += b;
 67110145263                this.m2[ind] += (r * r) + (g * g) + (b * b);
 264            }
 10265        }
 266
 267        /// <summary>
 268        /// Converts the histogram into moments so that we can rapidly calculate
 269        /// the sums of the above quantities over any desired box.
 270        /// </summary>
 271        private void Get3DMoments()
 272        {
 10273            Span<long> area = stackalloc long[IndexCount];
 10274            Span<long> areaR = stackalloc long[IndexCount];
 10275            Span<long> areaG = stackalloc long[IndexCount];
 10276            Span<long> areaB = stackalloc long[IndexCount];
 10277            Span<double> area2 = stackalloc double[IndexCount];
 278
 2580279            for (int r = 1; r < IndexCount; r++)
 280            {
 1280281                area.Clear();
 1280282                areaR.Clear();
 1280283                areaG.Clear();
 1280284                areaB.Clear();
 1280285                area2.Clear();
 286
 330240287                for (int g = 1; g < IndexCount; g++)
 288                {
 163840289                    long line = 0;
 163840290                    long lineR = 0;
 163840291                    long lineG = 0;
 163840292                    long lineB = 0;
 163840293                    double line2 = 0;
 294
 42270720295                    for (int b = 1; b < IndexCount; b++)
 296                    {
 20971520297                        int ind1 = GetIndex(r, g, b);
 298
 20971520299                        line += this.vwt[ind1];
 20971520300                        lineR += this.vmr[ind1];
 20971520301                        lineG += this.vmg[ind1];
 20971520302                        lineB += this.vmb[ind1];
 20971520303                        line2 += this.m2[ind1];
 304
 20971520305                        area[b] += line;
 20971520306                        areaR[b] += lineR;
 20971520307                        areaG[b] += lineG;
 20971520308                        areaB[b] += lineB;
 20971520309                        area2[b] += line2;
 310
 20971520311                        int ind2 = ind1 - GetIndex(1, 0, 0);
 312
 20971520313                        this.vwt[ind1] = this.vwt[ind2] + area[b];
 20971520314                        this.vmr[ind1] = this.vmr[ind2] + areaR[b];
 20971520315                        this.vmg[ind1] = this.vmg[ind2] + areaG[b];
 20971520316                        this.vmb[ind1] = this.vmb[ind2] + areaB[b];
 20971520317                        this.m2[ind1] = this.m2[ind2] + area2[b];
 318                    }
 319                }
 320            }
 10321        }
 322
 323        /// <summary>
 324        /// Computes the weighted variance of a box.
 325        /// </summary>
 326        /// <param name="cube">The cube.</param>
 327        /// <returns>The result.</returns>
 328        private double Variance(Box cube)
 329        {
 2541330            double dr = Volume(cube, this.vmr);
 2541331            double dg = Volume(cube, this.vmg);
 2541332            double db = Volume(cube, this.vmb);
 333
 2541334            double xx =
 2541335                this.m2[cube.Index111]
 2541336                - this.m2[cube.Index110]
 2541337                - this.m2[cube.Index101]
 2541338                + this.m2[cube.Index100]
 2541339                - this.m2[cube.Index011]
 2541340                + this.m2[cube.Index010]
 2541341                + this.m2[cube.Index001]
 2541342                - this.m2[cube.Index000];
 343
 2541344            return xx - (((dr * dr) + (dg * dg) + (db * db)) / Volume(cube, this.vwt));
 345        }
 346
 347        /// <summary>
 348        /// We want to minimize the sum of the variances of two sub-boxes.
 349        /// The sum(c^2) terms can be ignored since their sum over both sub-boxes
 350        /// is the same (the sum for the whole box) no matter where we split.
 351        /// The remaining terms have a minus sign in the variance formula,
 352        /// so we drop the minus sign and maximize the sum of the two terms.
 353        /// </summary>
 354        /// <param name="cube">The cube.</param>
 355        /// <param name="direction">The direction.</param>
 356        /// <param name="first">The first position.</param>
 357        /// <param name="last">The last position.</param>
 358        /// <param name="cut">The cutting point.</param>
 359        /// <param name="wholeR">The whole red.</param>
 360        /// <param name="wholeG">The whole green.</param>
 361        /// <param name="wholeB">The whole blue.</param>
 362        /// <param name="wholeW">The whole weight.</param>
 363        /// <returns>The result.</returns>
 364        private double Maximize(Box cube, int direction, int first, int last, out int cut, double wholeR, double wholeG,
 365        {
 6888366            long baseR = Bottom(cube, direction, this.vmr);
 6888367            long baseG = Bottom(cube, direction, this.vmg);
 6888368            long baseB = Bottom(cube, direction, this.vmb);
 6888369            long baseW = Bottom(cube, direction, this.vwt);
 370
 6888371            double max = 0.0;
 6888372            cut = -1;
 373
 1114880374            for (int i = first; i < last; i++)
 375            {
 550552376                double halfW1 = baseW + Top(cube, direction, i, this.vwt);
 377
 550552378                if (halfW1 == 0)
 379                {
 380                    continue;
 381                }
 382
 473296383                double halfW2 = wholeW - halfW1;
 384
 473296385                if (halfW2 == 0)
 386                {
 387                    continue;
 388                }
 389
 20923390                double halfR1 = baseR + Top(cube, direction, i, this.vmr);
 20923391                double halfG1 = baseG + Top(cube, direction, i, this.vmg);
 20923392                double halfB1 = baseB + Top(cube, direction, i, this.vmb);
 393
 20923394                double temp = ((halfR1 * halfR1) + (halfG1 * halfG1) + (halfB1 * halfB1)) / halfW1;
 395
 20923396                double halfR2 = wholeR - halfR1;
 20923397                double halfG2 = wholeG - halfG1;
 20923398                double halfB2 = wholeB - halfB1;
 399
 20923400                temp += ((halfR2 * halfR2) + (halfG2 * halfG2) + (halfB2 * halfB2)) / halfW2;
 401
 20923402                if (temp > max)
 403                {
 5888404                    max = temp;
 5888405                    cut = i;
 406                }
 407            }
 408
 6888409            return max;
 410        }
 411
 412        /// <summary>
 413        /// Cuts a box.
 414        /// </summary>
 415        /// <param name="set1">The first set.</param>
 416        /// <param name="set2">The second set.</param>
 417        /// <returns>Returns a value indicating whether the box has been split.</returns>
 418        private bool Cut(Box set1, Box set2)
 419        {
 2296420            double wholeR = Volume(set1, this.vmr);
 2296421            double wholeG = Volume(set1, this.vmg);
 2296422            double wholeB = Volume(set1, this.vmb);
 2296423            double wholeW = Volume(set1, this.vwt);
 424
 425            int cutr;
 426            int cutg;
 427            int cutb;
 428
 2296429            double maxr = this.Maximize(set1, 2, set1.R0 + 1, set1.R1, out cutr, wholeR, wholeG, wholeB, wholeW);
 2296430            double maxg = this.Maximize(set1, 1, set1.G0 + 1, set1.G1, out cutg, wholeR, wholeG, wholeB, wholeW);
 2296431            double maxb = this.Maximize(set1, 0, set1.B0 + 1, set1.B1, out cutb, wholeR, wholeG, wholeB, wholeW);
 432
 433            int dir;
 434
 2296435            if ((maxr >= maxg) && (maxr >= maxb))
 436            {
 1615437                dir = 2;
 438
 1615439                if (cutr < 0)
 440                {
 1025441                    return false;
 442                }
 443            }
 681444            else if ((maxg >= maxr) && (maxg >= maxb))
 445            {
 263446                dir = 1;
 447            }
 448            else
 449            {
 418450                dir = 0;
 451            }
 452
 1271453            set2.R1 = set1.R1;
 1271454            set2.G1 = set1.G1;
 1271455            set2.B1 = set1.B1;
 456
 457            switch (dir)
 458            {
 459                // Red
 460                case 2:
 590461                    set2.R0 = set1.R1 = cutr;
 590462                    set2.G0 = set1.G0;
 590463                    set2.B0 = set1.B0;
 590464                    break;
 465
 466                // Green
 467                case 1:
 263468                    set2.G0 = set1.G1 = cutg;
 263469                    set2.R0 = set1.R0;
 263470                    set2.B0 = set1.B0;
 263471                    break;
 472
 473                // Blue
 474                case 0:
 418475                    set2.B0 = set1.B1 = cutb;
 418476                    set2.R0 = set1.R0;
 418477                    set2.G0 = set1.G0;
 478                    break;
 479            }
 480
 1271481            set1.Update();
 1271482            set2.Update();
 483
 1271484            return true;
 485        }
 486
 487        /// <summary>
 488        /// Marks a color space tag.
 489        /// </summary>
 490        /// <param name="cube">The cube.</param>
 491        /// <param name="label">A label.</param>
 492        private void Mark(Box cube, byte label)
 493        {
 143106494            for (int r = cube.R0 + 1; r <= cube.R1; r++)
 495            {
 9020672496                for (int g = cube.G0 + 1; g <= cube.G1; g++)
 497                {
 50823168498                    for (int b = cube.B0 + 1; b <= cube.B1; b++)
 499                    {
 20971520500                        this.tag[GetIndex(r, g, b)] = label;
 501                    }
 502                }
 503            }
 1281504        }
 505
 506        /// <summary>
 507        /// Builds the cube.
 508        /// </summary>
 509        /// <param name="cube">The cube.</param>
 510        /// <param name="colorCount">The color count.</param>
 511        private void BuildCube(out Box[] cube, ref int colorCount)
 512        {
 10513            cube = new Box[colorCount];
 10514            double[] vv = new double[colorCount];
 515
 5140516            for (int i = 0; i < colorCount; i++)
 517            {
 2560518                cube[i] = new Box();
 519            }
 520
 10521            cube[0].R0 = cube[0].G0 = cube[0].B0 = 0;
 10522            cube[0].R1 = cube[0].G1 = cube[0].B1 = IndexCount - 1;
 10523            cube[0].Update();
 524
 10525            int next = 0;
 526
 4592527            for (int i = 1; i < colorCount; i++)
 528            {
 2296529                if (this.Cut(cube[next], cube[i]))
 530                {
 1271531                    vv[next] = cube[next].Volume > 1 ? this.Variance(cube[next]) : 0.0;
 1271532                    vv[i] = cube[i].Volume > 1 ? this.Variance(cube[i]) : 0.0;
 533                }
 534                else
 535                {
 1025536                    vv[next] = 0.0;
 1025537                    i--;
 538                }
 539
 2296540                next = 0;
 541
 2296542                double temp = vv[0];
 460016543                for (int k = 1; k <= i; k++)
 544                {
 227712545                    if (vv[k] > temp)
 546                    {
 2223547                        temp = vv[k];
 2223548                        next = k;
 549                    }
 550                }
 551
 2296552                if (temp <= 0.0)
 553                {
 10554                    colorCount = i + 1;
 10555                    break;
 556                }
 557            }
 0558        }
 559
 560        /// <summary>
 561        /// Generates the quantized result.
 562        /// </summary>
 563        /// <param name="image">The image.</param>
 564        /// <param name="colorCount">The color count.</param>
 565        /// <param name="cube">The cube.</param>
 566        /// <returns>The result.</returns>
 567        private ColorQuantizerResult GenerateResult(byte[] image, int colorCount, Box[] cube)
 568        {
 10569            var quantizedImage = new ColorQuantizerResult(image.Length / 4, colorCount);
 570
 2582571            for (int k = 0; k < colorCount; k++)
 572            {
 1281573                this.Mark(cube[k], (byte)k);
 574
 1281575                double weight = Volume(cube[k], this.vwt);
 576
 1281577                if (weight != 0)
 578                {
 1281579                    quantizedImage.Palette[(k * 4) + 3] = 0xff;
 1281580                    quantizedImage.Palette[(k * 4) + 2] = (byte)(Volume(cube[k], this.vmr) / weight);
 1281581                    quantizedImage.Palette[(k * 4) + 1] = (byte)(Volume(cube[k], this.vmg) / weight);
 1281582                    quantizedImage.Palette[k * 4] = (byte)(Volume(cube[k], this.vmb) / weight);
 583                }
 584                else
 585                {
 0586                    quantizedImage.Palette[(k * 4) + 3] = 0xff;
 0587                    quantizedImage.Palette[(k * 4) + 2] = 0;
 0588                    quantizedImage.Palette[(k * 4) + 1] = 0;
 0589                    quantizedImage.Palette[k * 4] = 0;
 590                }
 591            }
 592
 134220310593            for (int i = 0; i < image.Length / 4; i++)
 594            {
 67110145595                int r = image[(i * 4) + 2] >> (8 - IndexBits);
 67110145596                int g = image[(i * 4) + 1] >> (8 - IndexBits);
 67110145597                int b = image[i * 4] >> (8 - IndexBits);
 598
 67110145599                int ind = GetIndex(r + 1, g + 1, b + 1);
 600
 67110145601                quantizedImage.Bytes[i] = this.tag[ind];
 602            }
 603
 10604            return quantizedImage;
 605        }
 606    }
 607}
 608#endif
 609

https://raw.githubusercontent.com/JeremyAnsel/JeremyAnsel.ColorQuant/6d79217e72af9e3af1a8a29c606732adac1e8d87/JeremyAnsel.ColorQuant/JeremyAnsel.ColorQuant/WuColorQuantizer2.cs

#LineLine coverage
 1// <copyright file="WuColorQuantizer2.cs" company="Jérémy Ansel">
 2// Copyright (c) 2014-2019 Jérémy Ansel
 3// </copyright>
 4// <license>
 5// Licensed under the MIT license. See LICENSE.txt
 6// </license>
 7
 8#if NET8_0_OR_GREATER
 9namespace JeremyAnsel.ColorQuant
 10{
 11    using System;
 12    using System.Diagnostics.CodeAnalysis;
 13    using System.Runtime.Intrinsics;
 14
 15    /// <summary>
 16    /// A Wu's color quantizer.
 17    /// </summary>
 18    /// <remarks>
 19    /// <para>
 20    /// Based on C Implementation of Xiaolin Wu's Color Quantizer (v. 2)
 21    /// (see Graphics Gems volume II, pages 126-133)
 22    /// (<see href="http://www.ece.mcmaster.ca/~xwu/cq.c"/>).
 23    /// </para>
 24    /// <para>
 25    /// Algorithm: Greedy orthogonal bipartition of RGB space for variance
 26    /// minimization aided by inclusion-exclusion tricks.
 27    /// For speed no nearest neighbor search is done. Slightly
 28    /// better performance can be expected by more sophisticated
 29    /// but more expensive versions.
 30    /// </para>
 31    /// </remarks>
 32    [SuppressMessage("Microsoft.Naming", "CA1709:IdentifiersShouldBeCasedCorrectly", MessageId = "Wu", Justification = "
 33    public sealed class WuColorQuantizer : IColorQuantizer
 34    {
 35        /// <summary>
 36        /// The maximum color count for the quantize.
 37        /// </summary>
 38        public const int MaxColors = 256;
 39
 40        /// <summary>
 41        /// The index bits.
 42        /// </summary>
 43        private const int IndexBits = 7;
 44
 45        /// <summary>
 46        /// The index count.
 47        /// </summary>
 48        private const int IndexCount = (1 << IndexBits) + 1;
 49
 50        /// <summary>
 51        /// The table length.
 52        /// </summary>
 53        private const int TableLength = IndexCount * IndexCount * IndexCount;
 54
 155        private readonly Vector256<long>[] vm = new Vector256<long>[TableLength];
 56
 57        /// <summary>
 58        /// Moment of <c>c^2*P(c)</c>.
 59        /// </summary>
 160        private readonly double[] m2 = new double[TableLength];
 61
 62        /// <summary>
 63        /// Color space tag.
 64        /// </summary>
 165        private readonly byte[] tag = new byte[TableLength];
 66
 67        /// <summary>
 68        /// Quantizes an image.
 69        /// </summary>
 70        /// <param name="image">The image (XRGB).</param>
 71        /// <returns>The result.</returns>
 72        public ColorQuantizerResult Quantize(byte[] image)
 73        {
 1174            return this.Quantize(image, 256);
 75        }
 76
 77        /// <summary>
 78        /// Quantizes an image.
 79        /// </summary>
 80        /// <param name="image">The image (XRGB).</param>
 81        /// <param name="colorCount">The color count.</param>
 82        /// <returns>The result.</returns>
 83        public ColorQuantizerResult Quantize(byte[] image, int colorCount)
 84        {
 1485            if (image == null)
 86            {
 287                throw new ArgumentNullException(nameof(image));
 88            }
 89
 1290            if (colorCount < 1 || colorCount > 256)
 91            {
 292                throw new ArgumentOutOfRangeException(nameof(colorCount));
 93            }
 94
 1095            this.Clear();
 96
 1097            this.Build3DHistogram(image);
 1098            this.Get3DMoments();
 99
 100            Box[] cube;
 10101            this.BuildCube(out cube, ref colorCount);
 102
 10103            return this.GenerateResult(image, colorCount, cube);
 104        }
 105
 106        /// <summary>
 107        /// Gets an index.
 108        /// </summary>
 109        /// <param name="r">The red value.</param>
 110        /// <param name="g">The green value.</param>
 111        /// <param name="b">The blue value.</param>
 112        /// <returns>The index.</returns>
 113        private static int GetIndex(int r, int g, int b)
 114        {
 199337058115            return (r << (IndexBits * 2)) + (r << (IndexBits + 1)) + (g << IndexBits) + r + g + b;
 116        }
 117
 118        /// <summary>
 119        /// Computes sum over a box of any given statistic.
 120        /// </summary>
 121        /// <param name="cube">The cube.</param>
 122        /// <param name="moment">The moment.</param>
 123        /// <returns>The result.</returns>
 124        private static Vector256<long> Volume(Box cube, Vector256<long>[] moment)
 125        {
 6118126            return moment[cube.Index111]
 6118127               - moment[cube.Index110]
 6118128               - moment[cube.Index101]
 6118129               + moment[cube.Index100]
 6118130               - moment[cube.Index011]
 6118131               + moment[cube.Index010]
 6118132               + moment[cube.Index001]
 6118133               - moment[cube.Index000];
 134        }
 135
 136        /// <summary>
 137        /// Computes part of Volume(cube, moment) that doesn't depend on r1, g1, or b1 (depending on direction).
 138        /// </summary>
 139        /// <param name="cube">The cube.</param>
 140        /// <param name="direction">The direction.</param>
 141        /// <param name="moment">The moment.</param>
 142        /// <returns>The result.</returns>
 143        private static Vector256<long> Bottom(Box cube, int direction, Vector256<long>[] moment)
 144        {
 145            switch (direction)
 146            {
 147                // Red
 148                case 2:
 2296149                    return -moment[cube.Index011]
 2296150                        + moment[cube.Index010]
 2296151                        + moment[cube.Index001]
 2296152                        - moment[cube.Index000];
 153
 154                // Green
 155                case 1:
 2296156                    return -moment[cube.Index101]
 2296157                        + moment[cube.Index100]
 2296158                        + moment[cube.Index001]
 2296159                        - moment[cube.Index000];
 160
 161                // Blue
 162                case 0:
 2296163                    return -moment[cube.Index110]
 2296164                        + moment[cube.Index100]
 2296165                        + moment[cube.Index010]
 2296166                        - moment[cube.Index000];
 167
 168                default:
 0169                    throw new ArgumentOutOfRangeException(nameof(direction));
 170            }
 171        }
 172
 173        /// <summary>
 174        /// Computes remainder of Volume(cube, moment), substituting position for r1, g1, or b1 (depending on direction)
 175        /// </summary>
 176        /// <param name="cube">The cube.</param>
 177        /// <param name="direction">The direction.</param>
 178        /// <param name="position">The position.</param>
 179        /// <param name="moment">The moment.</param>
 180        /// <returns>The result.</returns>
 181        private static Vector256<long> Top(Box cube, int direction, int position, Vector256<long>[] moment)
 182        {
 183            switch (direction)
 184            {
 185                // Red
 186                case 2:
 139784187                    return moment[GetIndex(position, cube.G1, cube.B1)]
 139784188                       - moment[GetIndex(position, cube.G1, cube.B0)]
 139784189                       - moment[GetIndex(position, cube.G0, cube.B1)]
 139784190                       + moment[GetIndex(position, cube.G0, cube.B0)];
 191
 192                // Green
 193                case 1:
 204424194                    return moment[GetIndex(cube.R1, position, cube.B1)]
 204424195                       - moment[GetIndex(cube.R1, position, cube.B0)]
 204424196                       - moment[GetIndex(cube.R0, position, cube.B1)]
 204424197                       + moment[GetIndex(cube.R0, position, cube.B0)];
 198
 199                // Blue
 200                case 0:
 206344201                    return moment[GetIndex(cube.R1, cube.G1, position)]
 206344202                       - moment[GetIndex(cube.R1, cube.G0, position)]
 206344203                       - moment[GetIndex(cube.R0, cube.G1, position)]
 206344204                       + moment[GetIndex(cube.R0, cube.G0, position)];
 205
 206                default:
 0207                    throw new ArgumentOutOfRangeException(nameof(direction));
 208            }
 209        }
 210
 211        /// <summary>
 212        /// Clears the tables.
 213        /// </summary>
 214        private void Clear()
 215        {
 10216            this.vm.AsSpan().Clear();
 10217            this.m2.AsSpan().Clear();
 10218            this.tag.AsSpan().Clear();
 10219        }
 220
 221        /// <summary>
 222        /// Builds a 3-D color histogram of <c>counts, r/g/b, c^2</c>.
 223        /// </summary>
 224        /// <param name="image">The image.</param>
 225        private void Build3DHistogram(byte[] image)
 226        {
 134220310227            for (int i = 0; i < image.Length; i += 4)
 228            {
 67110145229                int r = image[i + 2];
 67110145230                int g = image[i + 1];
 67110145231                int b = image[i];
 232
 67110145233                int inr = r >> (8 - IndexBits);
 67110145234                int ing = g >> (8 - IndexBits);
 67110145235                int inb = b >> (8 - IndexBits);
 236
 67110145237                int ind = GetIndex(inr + 1, ing + 1, inb + 1);
 238
 67110145239                this.vm[ind] += Vector256.Create(r, g, b, 1);
 67110145240                this.m2[ind] += (r * r) + (g * g) + (b * b);
 241            }
 10242        }
 243
 244        /// <summary>
 245        /// Converts the histogram into moments so that we can rapidly calculate
 246        /// the sums of the above quantities over any desired box.
 247        /// </summary>
 248        private void Get3DMoments()
 249        {
 10250            Span<Vector256<long>> area = stackalloc Vector256<long>[IndexCount];
 10251            Span<double> area2 = stackalloc double[IndexCount];
 252
 2580253            for (int r = 1; r < IndexCount; r++)
 254            {
 1280255                area.Clear();
 1280256                area2.Clear();
 257
 330240258                for (int g = 1; g < IndexCount; g++)
 259                {
 163840260                    Vector256<long> line = Vector256<long>.Zero;
 163840261                    double line2 = 0;
 262
 42270720263                    for (int b = 1; b < IndexCount; b++)
 264                    {
 20971520265                        int ind1 = GetIndex(r, g, b);
 266
 20971520267                        line += this.vm[ind1];
 20971520268                        line2 += this.m2[ind1];
 269
 20971520270                        area[b] += line;
 20971520271                        area2[b] += line2;
 272
 20971520273                        int ind2 = ind1 - GetIndex(1, 0, 0);
 274
 20971520275                        this.vm[ind1] = this.vm[ind2] + area[b];
 20971520276                        this.m2[ind1] = this.m2[ind2] + area2[b];
 277                    }
 278                }
 279            }
 10280        }
 281
 282        /// <summary>
 283        /// Computes the weighted variance of a box.
 284        /// </summary>
 285        /// <param name="cube">The cube.</param>
 286        /// <returns>The result.</returns>
 287        private double Variance(Box cube)
 288        {
 2541289            Vector256<long> d = Volume(cube, this.vm);
 290
 2541291            double xx =
 2541292                this.m2[cube.Index111]
 2541293                - this.m2[cube.Index110]
 2541294                - this.m2[cube.Index101]
 2541295                + this.m2[cube.Index100]
 2541296                - this.m2[cube.Index011]
 2541297                + this.m2[cube.Index010]
 2541298                + this.m2[cube.Index001]
 2541299                - this.m2[cube.Index000];
 300
 2541301            Vector256<double> dv = Vector256.Create((double)d[0], (double)d[1], (double)d[2], 0.0);
 2541302            return xx - (Vector256.Dot(dv, dv) / d[3]);
 303        }
 304
 305        /// <summary>
 306        /// We want to minimize the sum of the variances of two sub-boxes.
 307        /// The sum(c^2) terms can be ignored since their sum over both sub-boxes
 308        /// is the same (the sum for the whole box) no matter where we split.
 309        /// The remaining terms have a minus sign in the variance formula,
 310        /// so we drop the minus sign and maximize the sum of the two terms.
 311        /// </summary>
 312        /// <param name="cube">The cube.</param>
 313        /// <param name="direction">The direction.</param>
 314        /// <param name="first">The first position.</param>
 315        /// <param name="last">The last position.</param>
 316        /// <param name="cut">The cutting point.</param>
 317        /// <param name="whole">The whole color and weight.</param>
 318        /// <returns>The result.</returns>
 319        private double Maximize(Box cube, int direction, int first, int last, out int cut, in Vector256<long> whole)
 320        {
 6888321            Vector256<long> bases = Bottom(cube, direction, this.vm);
 322
 6888323            double max = 0.0;
 6888324            cut = -1;
 325
 1114880326            for (int i = first; i < last; i++)
 327            {
 550552328                Vector256<long> halfs1 = bases + Top(cube, direction, i, this.vm);
 329
 550552330                if (halfs1[3] == 0)
 331                {
 332                    continue;
 333                }
 334
 473296335                Vector256<long> halfs2 = whole - halfs1;
 336
 473296337                if (halfs2[3] == 0)
 338                {
 339                    continue;
 340                }
 341
 20923342                Vector256<double> dv1 = Vector256.Create((double)halfs1[0], (double)halfs1[1], (double)halfs1[2], 0.0);
 20923343                double temp = Vector256.Dot(dv1, dv1) / halfs1[3];
 344
 20923345                Vector256<double> dv2 = Vector256.Create((double)halfs2[0], (double)halfs2[1], (double)halfs2[2], 0.0);
 20923346                temp += Vector256.Dot(dv2, dv2) / halfs2[3];
 347
 20923348                if (temp > max)
 349                {
 5888350                    max = temp;
 5888351                    cut = i;
 352                }
 353            }
 354
 6888355            return max;
 356        }
 357
 358        /// <summary>
 359        /// Cuts a box.
 360        /// </summary>
 361        /// <param name="set1">The first set.</param>
 362        /// <param name="set2">The second set.</param>
 363        /// <returns>Returns a value indicating whether the box has been split.</returns>
 364        private bool Cut(Box set1, Box set2)
 365        {
 2296366            Vector256<long> whole = Volume(set1, this.vm);
 367
 368            int cutr;
 369            int cutg;
 370            int cutb;
 371
 2296372            double maxr = this.Maximize(set1, 2, set1.R0 + 1, set1.R1, out cutr, whole);
 2296373            double maxg = this.Maximize(set1, 1, set1.G0 + 1, set1.G1, out cutg, whole);
 2296374            double maxb = this.Maximize(set1, 0, set1.B0 + 1, set1.B1, out cutb, whole);
 375
 376            int dir;
 377
 2296378            if ((maxr >= maxg) && (maxr >= maxb))
 379            {
 1615380                dir = 2;
 381
 1615382                if (cutr < 0)
 383                {
 1025384                    return false;
 385                }
 386            }
 681387            else if ((maxg >= maxr) && (maxg >= maxb))
 388            {
 263389                dir = 1;
 390            }
 391            else
 392            {
 418393                dir = 0;
 394            }
 395
 1271396            set2.R1 = set1.R1;
 1271397            set2.G1 = set1.G1;
 1271398            set2.B1 = set1.B1;
 399
 400            switch (dir)
 401            {
 402                // Red
 403                case 2:
 590404                    set2.R0 = set1.R1 = cutr;
 590405                    set2.G0 = set1.G0;
 590406                    set2.B0 = set1.B0;
 590407                    break;
 408
 409                // Green
 410                case 1:
 263411                    set2.G0 = set1.G1 = cutg;
 263412                    set2.R0 = set1.R0;
 263413                    set2.B0 = set1.B0;
 263414                    break;
 415
 416                // Blue
 417                case 0:
 418418                    set2.B0 = set1.B1 = cutb;
 418419                    set2.R0 = set1.R0;
 418420                    set2.G0 = set1.G0;
 421                    break;
 422            }
 423
 1271424            set1.Update();
 1271425            set2.Update();
 426
 1271427            return true;
 428        }
 429
 430        /// <summary>
 431        /// Marks a color space tag.
 432        /// </summary>
 433        /// <param name="cube">The cube.</param>
 434        /// <param name="label">A label.</param>
 435        private void Mark(Box cube, byte label)
 436        {
 143106437            for (int r = cube.R0 + 1; r <= cube.R1; r++)
 438            {
 9020672439                for (int g = cube.G0 + 1; g <= cube.G1; g++)
 440                {
 50823168441                    for (int b = cube.B0 + 1; b <= cube.B1; b++)
 442                    {
 20971520443                        this.tag[GetIndex(r, g, b)] = label;
 444                    }
 445                }
 446            }
 1281447        }
 448
 449        /// <summary>
 450        /// Builds the cube.
 451        /// </summary>
 452        /// <param name="cube">The cube.</param>
 453        /// <param name="colorCount">The color count.</param>
 454        private void BuildCube(out Box[] cube, ref int colorCount)
 455        {
 10456            cube = new Box[colorCount];
 10457            Span<double> vv = stackalloc double[colorCount];
 458
 5140459            for (int i = 0; i < colorCount; i++)
 460            {
 2560461                cube[i] = new Box();
 462            }
 463
 10464            cube[0].R0 = cube[0].G0 = cube[0].B0 = 0;
 10465            cube[0].R1 = cube[0].G1 = cube[0].B1 = IndexCount - 1;
 10466            cube[0].Update();
 467
 10468            int next = 0;
 469
 4592470            for (int i = 1; i < colorCount; i++)
 471            {
 2296472                if (this.Cut(cube[next], cube[i]))
 473                {
 1271474                    vv[next] = cube[next].Volume > 1 ? this.Variance(cube[next]) : 0.0;
 1271475                    vv[i] = cube[i].Volume > 1 ? this.Variance(cube[i]) : 0.0;
 476                }
 477                else
 478                {
 1025479                    vv[next] = 0.0;
 1025480                    i--;
 481                }
 482
 2296483                next = 0;
 484
 2296485                double temp = vv[0];
 460016486                for (int k = 1; k <= i; k++)
 487                {
 227712488                    if (vv[k] > temp)
 489                    {
 2223490                        temp = vv[k];
 2223491                        next = k;
 492                    }
 493                }
 494
 2296495                if (temp <= 0.0)
 496                {
 10497                    colorCount = i + 1;
 10498                    break;
 499                }
 500            }
 0501        }
 502
 503        /// <summary>
 504        /// Generates the quantized result.
 505        /// </summary>
 506        /// <param name="image">The image.</param>
 507        /// <param name="colorCount">The color count.</param>
 508        /// <param name="cube">The cube.</param>
 509        /// <returns>The result.</returns>
 510        private ColorQuantizerResult GenerateResult(byte[] image, int colorCount, Box[] cube)
 511        {
 10512            var quantizedImage = new ColorQuantizerResult(image.Length / 4, colorCount);
 513
 2582514            for (int k = 0; k < colorCount; k++)
 515            {
 1281516                this.Mark(cube[k], (byte)k);
 517
 1281518                Vector256<long> weight = Volume(cube[k], this.vm);
 1281519                double w = weight[3];
 520
 1281521                if (w != 0)
 522                {
 1281523                    quantizedImage.Palette[(k * 4) + 3] = 0xff;
 1281524                    quantizedImage.Palette[(k * 4) + 2] = (byte)(weight[0] / w);
 1281525                    quantizedImage.Palette[(k * 4) + 1] = (byte)(weight[1] / w);
 1281526                    quantizedImage.Palette[k * 4] = (byte)(weight[2] / w);
 527                }
 528                else
 529                {
 0530                    quantizedImage.Palette[(k * 4) + 3] = 0xff;
 0531                    quantizedImage.Palette[(k * 4) + 2] = 0;
 0532                    quantizedImage.Palette[(k * 4) + 1] = 0;
 0533                    quantizedImage.Palette[k * 4] = 0;
 534                }
 535            }
 536
 134220310537            for (int i = 0; i < image.Length / 4; i++)
 538            {
 67110145539                int r = image[(i * 4) + 2] >> (8 - IndexBits);
 67110145540                int g = image[(i * 4) + 1] >> (8 - IndexBits);
 67110145541                int b = image[i * 4] >> (8 - IndexBits);
 542
 67110145543                int ind = GetIndex(r + 1, g + 1, b + 1);
 544
 67110145545                quantizedImage.Bytes[i] = this.tag[ind];
 546            }
 547
 10548            return quantizedImage;
 549        }
 550    }
 551}
 552#endif
 553
Generated by: ReportGenerator 5.4.12.0
9/4/2025 - 6:55:30 PM
GitHub | reportgenerator.io

Methods/Properties

.ctor()
Quantize(System.Byte[])
Quantize(System.Byte[],System.Int32)
GetIndex(System.Int32,System.Int32,System.Int32)
Volume(JeremyAnsel.ColorQuant.Box,System.Int64[])
Bottom(JeremyAnsel.ColorQuant.Box,System.Int32,System.Int64[])
Top(JeremyAnsel.ColorQuant.Box,System.Int32,System.Int32,System.Int64[])
Clear()
Build3DHistogram(System.Byte[])
Get3DMoments()
Variance(JeremyAnsel.ColorQuant.Box)
Maximize(JeremyAnsel.ColorQuant.Box,System.Int32,System.Int32,System.Int32,System.Int32&,System.Double,System.Double,System.Double,System.Double)
Cut(JeremyAnsel.ColorQuant.Box,JeremyAnsel.ColorQuant.Box)
Mark(JeremyAnsel.ColorQuant.Box,System.Byte)
BuildCube(JeremyAnsel.ColorQuant.Box[]&,System.Int32&)
GenerateResult(System.Byte[],System.Int32,JeremyAnsel.ColorQuant.Box[])
.ctor()
Quantize(System.Byte[])
Quantize(System.Byte[],System.Int32)
GetIndex(System.Int32,System.Int32,System.Int32)
Volume(JeremyAnsel.ColorQuant.Box,System.Runtime.Intrinsics.Vector256`1<System.Int64>[])
Bottom(JeremyAnsel.ColorQuant.Box,System.Int32,System.Runtime.Intrinsics.Vector256`1<System.Int64>[])
Top(JeremyAnsel.ColorQuant.Box,System.Int32,System.Int32,System.Runtime.Intrinsics.Vector256`1<System.Int64>[])
Clear()
Build3DHistogram(System.Byte[])
Get3DMoments()
Variance(JeremyAnsel.ColorQuant.Box)
Maximize(JeremyAnsel.ColorQuant.Box,System.Int32,System.Int32,System.Int32,System.Int32&,System.Runtime.Intrinsics.Vector256`1<System.Int64>&)
Cut(JeremyAnsel.ColorQuant.Box,JeremyAnsel.ColorQuant.Box)
Mark(JeremyAnsel.ColorQuant.Box,System.Byte)
BuildCube(JeremyAnsel.ColorQuant.Box[]&,System.Int32&)
GenerateResult(System.Byte[],System.Int32,JeremyAnsel.ColorQuant.Box[])