mirror of https://github.com/SixLabors/ImageSharp
You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
795 lines
32 KiB
795 lines
32 KiB
// <copyright file="WuQuantizer.cs" company="James Jackson-South">
|
|
// Copyright © James Jackson-South and contributors.
|
|
// Licensed under the Apache License, Version 2.0.
|
|
// </copyright>
|
|
|
|
using System.Numerics;
|
|
|
|
namespace ImageProcessorCore.Quantizers
|
|
{
|
|
using System;
|
|
using System.Collections.Generic;
|
|
using System.Threading.Tasks;
|
|
|
|
/// <summary>
|
|
/// An implementation of Wu's color quantizer with alpha channel.
|
|
/// </summary>
|
|
/// <remarks>
|
|
/// <para>
|
|
/// Based on C Implementation of Xiaolin Wu's Color Quantizer (v. 2)
|
|
/// (see Graphics Gems volume II, pages 126-133)
|
|
/// (<see href="http://www.ece.mcmaster.ca/~xwu/cq.c"/>).
|
|
/// </para>
|
|
/// <para>
|
|
/// This adaptation is based on the excellent JeremyAnsel.ColorQuant by Jérémy Ansel
|
|
/// <see href="https://github.com/JeremyAnsel/JeremyAnsel.ColorQuant"/>
|
|
/// </para>
|
|
/// <para>
|
|
/// Algorithm: Greedy orthogonal bipartition of RGB space for variance
|
|
/// minimization aided by inclusion-exclusion tricks.
|
|
/// For speed no nearest neighbor search is done. Slightly
|
|
/// better performance can be expected by more sophisticated
|
|
/// but more expensive versions.
|
|
/// </para>
|
|
/// </remarks>
|
|
/// <typeparam name="TColor">The pixel format.</typeparam>
|
|
/// <typeparam name="TPacked">The packed format. <example>uint, long, float.</example></typeparam>
|
|
public sealed class WuQuantizer<TColor, TPacked> : IQuantizer<TColor, TPacked>
|
|
where TColor : IPackedVector<TPacked>
|
|
where TPacked : struct
|
|
{
|
|
/// <summary>
|
|
/// The epsilon for comparing floating point numbers.
|
|
/// </summary>
|
|
private const float Epsilon = 0.001f;
|
|
|
|
/// <summary>
|
|
/// The index bits.
|
|
/// </summary>
|
|
private const int IndexBits = 6;
|
|
|
|
/// <summary>
|
|
/// The index alpha bits.
|
|
/// </summary>
|
|
private const int IndexAlphaBits = 3;
|
|
|
|
/// <summary>
|
|
/// The index count.
|
|
/// </summary>
|
|
private const int IndexCount = (1 << IndexBits) + 1;
|
|
|
|
/// <summary>
|
|
/// The index alpha count.
|
|
/// </summary>
|
|
private const int IndexAlphaCount = (1 << IndexAlphaBits) + 1;
|
|
|
|
/// <summary>
|
|
/// The table length.
|
|
/// </summary>
|
|
private const int TableLength = IndexCount * IndexCount * IndexCount * IndexAlphaCount;
|
|
|
|
/// <summary>
|
|
/// Moment of <c>P(c)</c>.
|
|
/// </summary>
|
|
private readonly long[] vwt;
|
|
|
|
/// <summary>
|
|
/// Moment of <c>r*P(c)</c>.
|
|
/// </summary>
|
|
private readonly long[] vmr;
|
|
|
|
/// <summary>
|
|
/// Moment of <c>g*P(c)</c>.
|
|
/// </summary>
|
|
private readonly long[] vmg;
|
|
|
|
/// <summary>
|
|
/// Moment of <c>b*P(c)</c>.
|
|
/// </summary>
|
|
private readonly long[] vmb;
|
|
|
|
/// <summary>
|
|
/// Moment of <c>a*P(c)</c>.
|
|
/// </summary>
|
|
private readonly long[] vma;
|
|
|
|
/// <summary>
|
|
/// Moment of <c>c^2*P(c)</c>.
|
|
/// </summary>
|
|
private readonly double[] m2;
|
|
|
|
/// <summary>
|
|
/// Color space tag.
|
|
/// </summary>
|
|
private readonly byte[] tag;
|
|
|
|
/// <summary>
|
|
/// Initializes a new instance of the <see cref="WuQuantizer{T,TP}"/> class.
|
|
/// </summary>
|
|
public WuQuantizer()
|
|
{
|
|
this.vwt = new long[TableLength];
|
|
this.vmr = new long[TableLength];
|
|
this.vmg = new long[TableLength];
|
|
this.vmb = new long[TableLength];
|
|
this.vma = new long[TableLength];
|
|
this.m2 = new double[TableLength];
|
|
this.tag = new byte[TableLength];
|
|
}
|
|
|
|
/// <inheritdoc/>
|
|
public byte Threshold { get; set; }
|
|
|
|
/// <inheritdoc/>
|
|
public QuantizedImage<TColor, TPacked> Quantize(ImageBase<TColor, TPacked> image, int maxColors)
|
|
{
|
|
Guard.NotNull(image, nameof(image));
|
|
|
|
int colorCount = maxColors.Clamp(1, 256);
|
|
|
|
this.Clear();
|
|
|
|
using (PixelAccessor<TColor, TPacked> imagePixels = image.Lock())
|
|
{
|
|
this.Build3DHistogram(imagePixels);
|
|
this.Get3DMoments();
|
|
|
|
Box[] cube;
|
|
this.BuildCube(out cube, ref colorCount);
|
|
|
|
return this.GenerateResult(imagePixels, colorCount, cube);
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Gets an index.
|
|
/// </summary>
|
|
/// <param name="r">The red value.</param>
|
|
/// <param name="g">The green value.</param>
|
|
/// <param name="b">The blue value.</param>
|
|
/// <param name="a">The alpha value.</param>
|
|
/// <returns>The index.</returns>
|
|
private static int GetPaletteIndex(int r, int g, int b, int a)
|
|
{
|
|
return (r << ((IndexBits * 2) + IndexAlphaBits))
|
|
+ (r << (IndexBits + IndexAlphaBits + 1))
|
|
+ (g << (IndexBits + IndexAlphaBits))
|
|
+ (r << (IndexBits * 2))
|
|
+ (r << (IndexBits + 1))
|
|
+ (g << IndexBits)
|
|
+ ((r + g + b) << IndexAlphaBits)
|
|
+ r + g + b + a;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes sum over a box of any given statistic.
|
|
/// </summary>
|
|
/// <param name="cube">The cube.</param>
|
|
/// <param name="moment">The moment.</param>
|
|
/// <returns>The result.</returns>
|
|
private static double Volume(Box cube, long[] moment)
|
|
{
|
|
return moment[GetPaletteIndex(cube.R1, cube.G1, cube.B1, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G1, cube.B1, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G1, cube.B0, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G1, cube.B0, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G0, cube.B1, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G0, cube.B1, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G1, cube.B1, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G1, cube.B1, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A0)];
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes part of Volume(cube, moment) that doesn't depend on r1, g1, or b1 (depending on direction).
|
|
/// </summary>
|
|
/// <param name="cube">The cube.</param>
|
|
/// <param name="direction">The direction.</param>
|
|
/// <param name="moment">The moment.</param>
|
|
/// <returns>The result.</returns>
|
|
private static long Bottom(Box cube, int direction, long[] moment)
|
|
{
|
|
switch (direction)
|
|
{
|
|
// Red
|
|
case 0:
|
|
return -moment[GetPaletteIndex(cube.R0, cube.G1, cube.B1, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G1, cube.B1, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A0)];
|
|
|
|
// Green
|
|
case 1:
|
|
return -moment[GetPaletteIndex(cube.R1, cube.G0, cube.B1, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G0, cube.B1, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A0)];
|
|
|
|
// Blue
|
|
case 2:
|
|
return -moment[GetPaletteIndex(cube.R1, cube.G1, cube.B0, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G1, cube.B0, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A0)];
|
|
|
|
// Alpha
|
|
case 3:
|
|
return -moment[GetPaletteIndex(cube.R1, cube.G1, cube.B1, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G1, cube.B0, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G0, cube.B1, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G1, cube.B1, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A0)];
|
|
|
|
default:
|
|
throw new ArgumentOutOfRangeException(nameof(direction));
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes remainder of Volume(cube, moment), substituting position for r1, g1, or b1 (depending on direction).
|
|
/// </summary>
|
|
/// <param name="cube">The cube.</param>
|
|
/// <param name="direction">The direction.</param>
|
|
/// <param name="position">The position.</param>
|
|
/// <param name="moment">The moment.</param>
|
|
/// <returns>The result.</returns>
|
|
private static long Top(Box cube, int direction, int position, long[] moment)
|
|
{
|
|
switch (direction)
|
|
{
|
|
// Red
|
|
case 0:
|
|
return moment[GetPaletteIndex(position, cube.G1, cube.B1, cube.A1)]
|
|
- moment[GetPaletteIndex(position, cube.G1, cube.B1, cube.A0)]
|
|
- moment[GetPaletteIndex(position, cube.G1, cube.B0, cube.A1)]
|
|
+ moment[GetPaletteIndex(position, cube.G1, cube.B0, cube.A0)]
|
|
- moment[GetPaletteIndex(position, cube.G0, cube.B1, cube.A1)]
|
|
+ moment[GetPaletteIndex(position, cube.G0, cube.B1, cube.A0)]
|
|
+ moment[GetPaletteIndex(position, cube.G0, cube.B0, cube.A1)]
|
|
- moment[GetPaletteIndex(position, cube.G0, cube.B0, cube.A0)];
|
|
|
|
// Green
|
|
case 1:
|
|
return moment[GetPaletteIndex(cube.R1, position, cube.B1, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R1, position, cube.B1, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R1, position, cube.B0, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R1, position, cube.B0, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, position, cube.B1, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R0, position, cube.B1, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, position, cube.B0, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R0, position, cube.B0, cube.A0)];
|
|
|
|
// Blue
|
|
case 2:
|
|
return moment[GetPaletteIndex(cube.R1, cube.G1, position, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G1, position, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G0, position, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G0, position, cube.A0)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G1, position, cube.A1)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G1, position, cube.A0)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, position, cube.A1)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, position, cube.A0)];
|
|
|
|
// Alpha
|
|
case 3:
|
|
return moment[GetPaletteIndex(cube.R1, cube.G1, cube.B1, position)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G1, cube.B0, position)]
|
|
- moment[GetPaletteIndex(cube.R1, cube.G0, cube.B1, position)]
|
|
+ moment[GetPaletteIndex(cube.R1, cube.G0, cube.B0, position)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G1, cube.B1, position)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G1, cube.B0, position)]
|
|
+ moment[GetPaletteIndex(cube.R0, cube.G0, cube.B1, position)]
|
|
- moment[GetPaletteIndex(cube.R0, cube.G0, cube.B0, position)];
|
|
|
|
default:
|
|
throw new ArgumentOutOfRangeException(nameof(direction));
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Clears the tables.
|
|
/// </summary>
|
|
private void Clear()
|
|
{
|
|
Array.Clear(this.vwt, 0, TableLength);
|
|
Array.Clear(this.vmr, 0, TableLength);
|
|
Array.Clear(this.vmg, 0, TableLength);
|
|
Array.Clear(this.vmb, 0, TableLength);
|
|
Array.Clear(this.vma, 0, TableLength);
|
|
Array.Clear(this.m2, 0, TableLength);
|
|
|
|
Array.Clear(this.tag, 0, TableLength);
|
|
}
|
|
|
|
/// <summary>
|
|
/// Builds a 3-D color histogram of <c>counts, r/g/b, c^2</c>.
|
|
/// </summary>
|
|
/// <param name="pixels">The pixel accessor.</param>
|
|
private void Build3DHistogram(PixelAccessor<TColor, TPacked> pixels)
|
|
{
|
|
for (int y = 0; y < pixels.Height; y++)
|
|
{
|
|
for (int x = 0; x < pixels.Width; x++)
|
|
{
|
|
// Colors are expected in r->g->b->a format
|
|
Color color = new Color(pixels[x, y].ToVector4());
|
|
|
|
byte r = color.R;
|
|
byte g = color.G;
|
|
byte b = color.B;
|
|
byte a = color.A;
|
|
|
|
int inr = r >> (8 - IndexBits);
|
|
int ing = g >> (8 - IndexBits);
|
|
int inb = b >> (8 - IndexBits);
|
|
int ina = a >> (8 - IndexAlphaBits);
|
|
|
|
int ind = GetPaletteIndex(inr + 1, ing + 1, inb + 1, ina + 1);
|
|
|
|
this.vwt[ind]++;
|
|
this.vmr[ind] += r;
|
|
this.vmg[ind] += g;
|
|
this.vmb[ind] += b;
|
|
this.vma[ind] += a;
|
|
this.m2[ind] += (r * r) + (g * g) + (b * b) + (a * a);
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Converts the histogram into moments so that we can rapidly calculate
|
|
/// the sums of the above quantities over any desired box.
|
|
/// </summary>
|
|
private void Get3DMoments()
|
|
{
|
|
long[] volume = new long[IndexCount * IndexAlphaCount];
|
|
long[] volumeR = new long[IndexCount * IndexAlphaCount];
|
|
long[] volumeG = new long[IndexCount * IndexAlphaCount];
|
|
long[] volumeB = new long[IndexCount * IndexAlphaCount];
|
|
long[] volumeA = new long[IndexCount * IndexAlphaCount];
|
|
double[] volume2 = new double[IndexCount * IndexAlphaCount];
|
|
|
|
long[] area = new long[IndexAlphaCount];
|
|
long[] areaR = new long[IndexAlphaCount];
|
|
long[] areaG = new long[IndexAlphaCount];
|
|
long[] areaB = new long[IndexAlphaCount];
|
|
long[] areaA = new long[IndexAlphaCount];
|
|
double[] area2 = new double[IndexAlphaCount];
|
|
|
|
for (int r = 1; r < IndexCount; r++)
|
|
{
|
|
Array.Clear(volume, 0, IndexCount * IndexAlphaCount);
|
|
Array.Clear(volumeR, 0, IndexCount * IndexAlphaCount);
|
|
Array.Clear(volumeG, 0, IndexCount * IndexAlphaCount);
|
|
Array.Clear(volumeB, 0, IndexCount * IndexAlphaCount);
|
|
Array.Clear(volumeA, 0, IndexCount * IndexAlphaCount);
|
|
Array.Clear(volume2, 0, IndexCount * IndexAlphaCount);
|
|
|
|
for (int g = 1; g < IndexCount; g++)
|
|
{
|
|
Array.Clear(area, 0, IndexAlphaCount);
|
|
Array.Clear(areaR, 0, IndexAlphaCount);
|
|
Array.Clear(areaG, 0, IndexAlphaCount);
|
|
Array.Clear(areaB, 0, IndexAlphaCount);
|
|
Array.Clear(areaA, 0, IndexAlphaCount);
|
|
Array.Clear(area2, 0, IndexAlphaCount);
|
|
|
|
for (int b = 1; b < IndexCount; b++)
|
|
{
|
|
long line = 0;
|
|
long lineR = 0;
|
|
long lineG = 0;
|
|
long lineB = 0;
|
|
long lineA = 0;
|
|
double line2 = 0;
|
|
|
|
for (int a = 1; a < IndexAlphaCount; a++)
|
|
{
|
|
int ind1 = GetPaletteIndex(r, g, b, a);
|
|
|
|
line += this.vwt[ind1];
|
|
lineR += this.vmr[ind1];
|
|
lineG += this.vmg[ind1];
|
|
lineB += this.vmb[ind1];
|
|
lineA += this.vma[ind1];
|
|
line2 += this.m2[ind1];
|
|
|
|
area[a] += line;
|
|
areaR[a] += lineR;
|
|
areaG[a] += lineG;
|
|
areaB[a] += lineB;
|
|
areaA[a] += lineA;
|
|
area2[a] += line2;
|
|
|
|
int inv = (b * IndexAlphaCount) + a;
|
|
|
|
volume[inv] += area[a];
|
|
volumeR[inv] += areaR[a];
|
|
volumeG[inv] += areaG[a];
|
|
volumeB[inv] += areaB[a];
|
|
volumeA[inv] += areaA[a];
|
|
volume2[inv] += area2[a];
|
|
|
|
int ind2 = ind1 - GetPaletteIndex(1, 0, 0, 0);
|
|
|
|
this.vwt[ind1] = this.vwt[ind2] + volume[inv];
|
|
this.vmr[ind1] = this.vmr[ind2] + volumeR[inv];
|
|
this.vmg[ind1] = this.vmg[ind2] + volumeG[inv];
|
|
this.vmb[ind1] = this.vmb[ind2] + volumeB[inv];
|
|
this.vma[ind1] = this.vma[ind2] + volumeA[inv];
|
|
this.m2[ind1] = this.m2[ind2] + volume2[inv];
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Computes the weighted variance of a box cube.
|
|
/// </summary>
|
|
/// <param name="cube">The cube.</param>
|
|
/// <returns>The <see cref="double"/>.</returns>
|
|
private double Variance(Box cube)
|
|
{
|
|
double dr = Volume(cube, this.vmr);
|
|
double dg = Volume(cube, this.vmg);
|
|
double db = Volume(cube, this.vmb);
|
|
double da = Volume(cube, this.vma);
|
|
|
|
double xx =
|
|
this.m2[GetPaletteIndex(cube.R1, cube.G1, cube.B1, cube.A1)]
|
|
- this.m2[GetPaletteIndex(cube.R1, cube.G1, cube.B1, cube.A0)]
|
|
- this.m2[GetPaletteIndex(cube.R1, cube.G1, cube.B0, cube.A1)]
|
|
+ this.m2[GetPaletteIndex(cube.R1, cube.G1, cube.B0, cube.A0)]
|
|
- this.m2[GetPaletteIndex(cube.R1, cube.G0, cube.B1, cube.A1)]
|
|
+ this.m2[GetPaletteIndex(cube.R1, cube.G0, cube.B1, cube.A0)]
|
|
+ this.m2[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A1)]
|
|
- this.m2[GetPaletteIndex(cube.R1, cube.G0, cube.B0, cube.A0)]
|
|
- this.m2[GetPaletteIndex(cube.R0, cube.G1, cube.B1, cube.A1)]
|
|
+ this.m2[GetPaletteIndex(cube.R0, cube.G1, cube.B1, cube.A0)]
|
|
+ this.m2[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A1)]
|
|
- this.m2[GetPaletteIndex(cube.R0, cube.G1, cube.B0, cube.A0)]
|
|
+ this.m2[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A1)]
|
|
- this.m2[GetPaletteIndex(cube.R0, cube.G0, cube.B1, cube.A0)]
|
|
- this.m2[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A1)]
|
|
+ this.m2[GetPaletteIndex(cube.R0, cube.G0, cube.B0, cube.A0)];
|
|
|
|
return xx - (((dr * dr) + (dg * dg) + (db * db) + (da * da)) / Volume(cube, this.vwt));
|
|
}
|
|
|
|
/// <summary>
|
|
/// We want to minimize the sum of the variances of two sub-boxes.
|
|
/// The sum(c^2) terms can be ignored since their sum over both sub-boxes
|
|
/// is the same (the sum for the whole box) no matter where we split.
|
|
/// The remaining terms have a minus sign in the variance formula,
|
|
/// so we drop the minus sign and maximize the sum of the two terms.
|
|
/// </summary>
|
|
/// <param name="cube">The cube.</param>
|
|
/// <param name="direction">The direction.</param>
|
|
/// <param name="first">The first position.</param>
|
|
/// <param name="last">The last position.</param>
|
|
/// <param name="cut">The cutting point.</param>
|
|
/// <param name="wholeR">The whole red.</param>
|
|
/// <param name="wholeG">The whole green.</param>
|
|
/// <param name="wholeB">The whole blue.</param>
|
|
/// <param name="wholeA">The whole alpha.</param>
|
|
/// <param name="wholeW">The whole weight.</param>
|
|
/// <returns>The <see cref="double"/>.</returns>
|
|
private double Maximize(Box cube, int direction, int first, int last, out int cut, double wholeR, double wholeG, double wholeB, double wholeA, double wholeW)
|
|
{
|
|
long baseR = Bottom(cube, direction, this.vmr);
|
|
long baseG = Bottom(cube, direction, this.vmg);
|
|
long baseB = Bottom(cube, direction, this.vmb);
|
|
long baseA = Bottom(cube, direction, this.vma);
|
|
long baseW = Bottom(cube, direction, this.vwt);
|
|
|
|
double max = 0.0;
|
|
cut = -1;
|
|
|
|
for (int i = first; i < last; i++)
|
|
{
|
|
double halfR = baseR + Top(cube, direction, i, this.vmr);
|
|
double halfG = baseG + Top(cube, direction, i, this.vmg);
|
|
double halfB = baseB + Top(cube, direction, i, this.vmb);
|
|
double halfA = baseA + Top(cube, direction, i, this.vma);
|
|
double halfW = baseW + Top(cube, direction, i, this.vwt);
|
|
|
|
double temp;
|
|
|
|
if (Math.Abs(halfW) < Epsilon)
|
|
{
|
|
continue;
|
|
}
|
|
|
|
temp = ((halfR * halfR) + (halfG * halfG) + (halfB * halfB) + (halfA * halfA)) / halfW;
|
|
|
|
halfR = wholeR - halfR;
|
|
halfG = wholeG - halfG;
|
|
halfB = wholeB - halfB;
|
|
halfA = wholeA - halfA;
|
|
halfW = wholeW - halfW;
|
|
|
|
if (Math.Abs(halfW) < Epsilon)
|
|
{
|
|
continue;
|
|
}
|
|
|
|
temp += ((halfR * halfR) + (halfG * halfG) + (halfB * halfB) + (halfA * halfA)) / halfW;
|
|
|
|
if (temp > max)
|
|
{
|
|
max = temp;
|
|
cut = i;
|
|
}
|
|
}
|
|
|
|
return max;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Cuts a box.
|
|
/// </summary>
|
|
/// <param name="set1">The first set.</param>
|
|
/// <param name="set2">The second set.</param>
|
|
/// <returns>Returns a value indicating whether the box has been split.</returns>
|
|
private bool Cut(Box set1, Box set2)
|
|
{
|
|
double wholeR = Volume(set1, this.vmr);
|
|
double wholeG = Volume(set1, this.vmg);
|
|
double wholeB = Volume(set1, this.vmb);
|
|
double wholeA = Volume(set1, this.vma);
|
|
double wholeW = Volume(set1, this.vwt);
|
|
|
|
int cutr;
|
|
int cutg;
|
|
int cutb;
|
|
int cuta;
|
|
|
|
double maxr = this.Maximize(set1, 0, set1.R0 + 1, set1.R1, out cutr, wholeR, wholeG, wholeB, wholeA, wholeW);
|
|
double maxg = this.Maximize(set1, 1, set1.G0 + 1, set1.G1, out cutg, wholeR, wholeG, wholeB, wholeA, wholeW);
|
|
double maxb = this.Maximize(set1, 2, set1.B0 + 1, set1.B1, out cutb, wholeR, wholeG, wholeB, wholeA, wholeW);
|
|
double maxa = this.Maximize(set1, 3, set1.A0 + 1, set1.A1, out cuta, wholeR, wholeG, wholeB, wholeA, wholeW);
|
|
|
|
int dir;
|
|
|
|
if ((maxr >= maxg) && (maxr >= maxb) && (maxr >= maxa))
|
|
{
|
|
dir = 0;
|
|
|
|
if (cutr < 0)
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
else if ((maxg >= maxr) && (maxg >= maxb) && (maxg >= maxa))
|
|
{
|
|
dir = 1;
|
|
}
|
|
else if ((maxb >= maxr) && (maxb >= maxg) && (maxb >= maxa))
|
|
{
|
|
dir = 2;
|
|
}
|
|
else
|
|
{
|
|
dir = 3;
|
|
}
|
|
|
|
set2.R1 = set1.R1;
|
|
set2.G1 = set1.G1;
|
|
set2.B1 = set1.B1;
|
|
set2.A1 = set1.A1;
|
|
|
|
switch (dir)
|
|
{
|
|
// Red
|
|
case 0:
|
|
set2.R0 = set1.R1 = cutr;
|
|
set2.G0 = set1.G0;
|
|
set2.B0 = set1.B0;
|
|
set2.A0 = set1.A0;
|
|
break;
|
|
|
|
// Green
|
|
case 1:
|
|
set2.G0 = set1.G1 = cutg;
|
|
set2.R0 = set1.R0;
|
|
set2.B0 = set1.B0;
|
|
set2.A0 = set1.A0;
|
|
break;
|
|
|
|
// Blue
|
|
case 2:
|
|
set2.B0 = set1.B1 = cutb;
|
|
set2.R0 = set1.R0;
|
|
set2.G0 = set1.G0;
|
|
set2.A0 = set1.A0;
|
|
break;
|
|
|
|
// Alpha
|
|
case 3:
|
|
set2.A0 = set1.A1 = cuta;
|
|
set2.R0 = set1.R0;
|
|
set2.G0 = set1.G0;
|
|
set2.B0 = set1.B0;
|
|
break;
|
|
}
|
|
|
|
set1.Volume = (set1.R1 - set1.R0) * (set1.G1 - set1.G0) * (set1.B1 - set1.B0) * (set1.A1 - set1.A0);
|
|
set2.Volume = (set2.R1 - set2.R0) * (set2.G1 - set2.G0) * (set2.B1 - set2.B0) * (set2.A1 - set2.A0);
|
|
|
|
return true;
|
|
}
|
|
|
|
/// <summary>
|
|
/// Marks a color space tag.
|
|
/// </summary>
|
|
/// <param name="cube">The cube.</param>
|
|
/// <param name="label">A label.</param>
|
|
private void Mark(Box cube, byte label)
|
|
{
|
|
for (int r = cube.R0 + 1; r <= cube.R1; r++)
|
|
{
|
|
for (int g = cube.G0 + 1; g <= cube.G1; g++)
|
|
{
|
|
for (int b = cube.B0 + 1; b <= cube.B1; b++)
|
|
{
|
|
for (int a = cube.A0 + 1; a <= cube.A1; a++)
|
|
{
|
|
this.tag[GetPaletteIndex(r, g, b, a)] = label;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Builds the cube.
|
|
/// </summary>
|
|
/// <param name="cube">The cube.</param>
|
|
/// <param name="colorCount">The color count.</param>
|
|
private void BuildCube(out Box[] cube, ref int colorCount)
|
|
{
|
|
cube = new Box[colorCount];
|
|
double[] vv = new double[colorCount];
|
|
|
|
for (int i = 0; i < colorCount; i++)
|
|
{
|
|
cube[i] = new Box();
|
|
}
|
|
|
|
cube[0].R0 = cube[0].G0 = cube[0].B0 = cube[0].A0 = 0;
|
|
cube[0].R1 = cube[0].G1 = cube[0].B1 = IndexCount - 1;
|
|
cube[0].A1 = IndexAlphaCount - 1;
|
|
|
|
int next = 0;
|
|
|
|
for (int i = 1; i < colorCount; i++)
|
|
{
|
|
if (this.Cut(cube[next], cube[i]))
|
|
{
|
|
vv[next] = cube[next].Volume > 1 ? this.Variance(cube[next]) : 0.0;
|
|
vv[i] = cube[i].Volume > 1 ? this.Variance(cube[i]) : 0.0;
|
|
}
|
|
else
|
|
{
|
|
vv[next] = 0.0;
|
|
i--;
|
|
}
|
|
|
|
next = 0;
|
|
|
|
double temp = vv[0];
|
|
for (int k = 1; k <= i; k++)
|
|
{
|
|
if (vv[k] > temp)
|
|
{
|
|
temp = vv[k];
|
|
next = k;
|
|
}
|
|
}
|
|
|
|
if (temp <= 0.0)
|
|
{
|
|
colorCount = i + 1;
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
/// <summary>
|
|
/// Generates the quantized result.
|
|
/// </summary>
|
|
/// <param name="imagePixels">The image pixels.</param>
|
|
/// <param name="colorCount">The color count.</param>
|
|
/// <param name="cube">The cube.</param>
|
|
/// <returns>The result.</returns>
|
|
private QuantizedImage<TColor, TPacked> GenerateResult(PixelAccessor<TColor, TPacked> imagePixels, int colorCount, Box[] cube)
|
|
{
|
|
List<TColor> pallette = new List<TColor>();
|
|
byte[] pixels = new byte[imagePixels.Width * imagePixels.Height];
|
|
int transparentIndex = -1;
|
|
int width = imagePixels.Width;
|
|
int height = imagePixels.Height;
|
|
|
|
for (int k = 0; k < colorCount; k++)
|
|
{
|
|
this.Mark(cube[k], (byte)k);
|
|
|
|
double weight = Volume(cube[k], this.vwt);
|
|
|
|
if (Math.Abs(weight) > Epsilon)
|
|
{
|
|
byte r = (byte)(Volume(cube[k], this.vmr) / weight);
|
|
byte g = (byte)(Volume(cube[k], this.vmg) / weight);
|
|
byte b = (byte)(Volume(cube[k], this.vmb) / weight);
|
|
byte a = (byte)(Volume(cube[k], this.vma) / weight);
|
|
|
|
TColor color = default(TColor);
|
|
color.PackFromVector4(new Vector4(r, g, b, a) / 255F);
|
|
|
|
if (color.Equals(default(TColor)))
|
|
{
|
|
transparentIndex = k;
|
|
}
|
|
|
|
pallette.Add(color);
|
|
}
|
|
else
|
|
{
|
|
pallette.Add(default(TColor));
|
|
transparentIndex = k;
|
|
}
|
|
}
|
|
|
|
Parallel.For(
|
|
0,
|
|
height,
|
|
Bootstrapper.Instance.ParallelOptions,
|
|
y =>
|
|
{
|
|
for (int x = 0; x < width; x++)
|
|
{
|
|
// Expected order r->g->b->a
|
|
Color color = new Color(imagePixels[x, y].ToVector4());
|
|
int r = color.R >> (8 - IndexBits);
|
|
int g = color.G >> (8 - IndexBits);
|
|
int b = color.B >> (8 - IndexBits);
|
|
int a = color.A >> (8 - IndexAlphaBits);
|
|
|
|
if (transparentIndex > -1 && color.A <= this.Threshold)
|
|
{
|
|
pixels[(y * width) + x] = (byte)transparentIndex;
|
|
continue;
|
|
}
|
|
|
|
int ind = GetPaletteIndex(r + 1, g + 1, b + 1, a + 1);
|
|
pixels[(y * width) + x] = this.tag[ind];
|
|
}
|
|
});
|
|
|
|
return new QuantizedImage<TColor, TPacked>(width, height, pallette.ToArray(), pixels, transparentIndex);
|
|
}
|
|
}
|
|
}
|