1 /* 2 * Copyright 2021 Google LLC 3 * 4 * Licensed under the Apache License, Version 2.0 (the "License"); 5 * you may not use this file except in compliance with the License. 6 * You may obtain a copy of the License at 7 * 8 * http://www.apache.org/licenses/LICENSE-2.0 9 * 10 * Unless required by applicable law or agreed to in writing, software 11 * distributed under the License is distributed on an "AS IS" BASIS, 12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 13 * See the License for the specific language governing permissions and 14 * limitations under the License. 15 */ 16 17 package com.google.ux.material.libmonet.quantize; 18 19 import com.google.ux.material.libmonet.utils.ColorUtils; 20 import java.util.ArrayList; 21 import java.util.LinkedHashMap; 22 import java.util.List; 23 import java.util.Map; 24 25 /** 26 * An image quantizer that divides the image's pixels into clusters by recursively cutting an RGB 27 * cube, based on the weight of pixels in each area of the cube. 28 * 29 * <p>The algorithm was described by Xiaolin Wu in Graphic Gems II, published in 1991. 30 */ 31 public final class QuantizerWu implements Quantizer { 32 int[] weights; 33 int[] momentsR; 34 int[] momentsG; 35 int[] momentsB; 36 double[] moments; 37 Box[] cubes; 38 39 // A histogram of all the input colors is constructed. It has the shape of a 40 // cube. The cube would be too large if it contained all 16 million colors: 41 // historical best practice is to use 5 bits of the 8 in each channel, 42 // reducing the histogram to a volume of ~32,000. 43 private static final int INDEX_BITS = 5; 44 private static final int INDEX_COUNT = 33; // ((1 << INDEX_BITS) + 1) 45 private static final int TOTAL_SIZE = 35937; // INDEX_COUNT * INDEX_COUNT * INDEX_COUNT 46 47 @Override quantize(int[] pixels, int colorCount)48 public QuantizerResult quantize(int[] pixels, int colorCount) { 49 QuantizerResult mapResult = new QuantizerMap().quantize(pixels, colorCount); 50 constructHistogram(mapResult.colorToCount); 51 createMoments(); 52 CreateBoxesResult createBoxesResult = createBoxes(colorCount); 53 List<Integer> colors = createResult(createBoxesResult.resultCount); 54 Map<Integer, Integer> resultMap = new LinkedHashMap<>(); 55 for (int color : colors) { 56 resultMap.put(color, 0); 57 } 58 return new QuantizerResult(resultMap); 59 } 60 getIndex(int r, int g, int b)61 static int getIndex(int r, int g, int b) { 62 return (r << (INDEX_BITS * 2)) + (r << (INDEX_BITS + 1)) + r + (g << INDEX_BITS) + g + b; 63 } 64 constructHistogram(Map<Integer, Integer> pixels)65 void constructHistogram(Map<Integer, Integer> pixels) { 66 weights = new int[TOTAL_SIZE]; 67 momentsR = new int[TOTAL_SIZE]; 68 momentsG = new int[TOTAL_SIZE]; 69 momentsB = new int[TOTAL_SIZE]; 70 moments = new double[TOTAL_SIZE]; 71 72 for (Map.Entry<Integer, Integer> pair : pixels.entrySet()) { 73 int pixel = pair.getKey(); 74 int count = pair.getValue(); 75 int red = ColorUtils.redFromArgb(pixel); 76 int green = ColorUtils.greenFromArgb(pixel); 77 int blue = ColorUtils.blueFromArgb(pixel); 78 int bitsToRemove = 8 - INDEX_BITS; 79 int iR = (red >> bitsToRemove) + 1; 80 int iG = (green >> bitsToRemove) + 1; 81 int iB = (blue >> bitsToRemove) + 1; 82 int index = getIndex(iR, iG, iB); 83 weights[index] += count; 84 momentsR[index] += (red * count); 85 momentsG[index] += (green * count); 86 momentsB[index] += (blue * count); 87 moments[index] += (count * ((red * red) + (green * green) + (blue * blue))); 88 } 89 } 90 createMoments()91 void createMoments() { 92 for (int r = 1; r < INDEX_COUNT; ++r) { 93 int[] area = new int[INDEX_COUNT]; 94 int[] areaR = new int[INDEX_COUNT]; 95 int[] areaG = new int[INDEX_COUNT]; 96 int[] areaB = new int[INDEX_COUNT]; 97 double[] area2 = new double[INDEX_COUNT]; 98 99 for (int g = 1; g < INDEX_COUNT; ++g) { 100 int line = 0; 101 int lineR = 0; 102 int lineG = 0; 103 int lineB = 0; 104 double line2 = 0.0; 105 for (int b = 1; b < INDEX_COUNT; ++b) { 106 int index = getIndex(r, g, b); 107 line += weights[index]; 108 lineR += momentsR[index]; 109 lineG += momentsG[index]; 110 lineB += momentsB[index]; 111 line2 += moments[index]; 112 113 area[b] += line; 114 areaR[b] += lineR; 115 areaG[b] += lineG; 116 areaB[b] += lineB; 117 area2[b] += line2; 118 119 int previousIndex = getIndex(r - 1, g, b); 120 weights[index] = weights[previousIndex] + area[b]; 121 momentsR[index] = momentsR[previousIndex] + areaR[b]; 122 momentsG[index] = momentsG[previousIndex] + areaG[b]; 123 momentsB[index] = momentsB[previousIndex] + areaB[b]; 124 moments[index] = moments[previousIndex] + area2[b]; 125 } 126 } 127 } 128 } 129 createBoxes(int maxColorCount)130 CreateBoxesResult createBoxes(int maxColorCount) { 131 cubes = new Box[maxColorCount]; 132 for (int i = 0; i < maxColorCount; i++) { 133 cubes[i] = new Box(); 134 } 135 double[] volumeVariance = new double[maxColorCount]; 136 Box firstBox = cubes[0]; 137 firstBox.r1 = INDEX_COUNT - 1; 138 firstBox.g1 = INDEX_COUNT - 1; 139 firstBox.b1 = INDEX_COUNT - 1; 140 141 int generatedColorCount = maxColorCount; 142 int next = 0; 143 for (int i = 1; i < maxColorCount; i++) { 144 if (cut(cubes[next], cubes[i])) { 145 volumeVariance[next] = (cubes[next].vol > 1) ? variance(cubes[next]) : 0.0; 146 volumeVariance[i] = (cubes[i].vol > 1) ? variance(cubes[i]) : 0.0; 147 } else { 148 volumeVariance[next] = 0.0; 149 i--; 150 } 151 152 next = 0; 153 154 double temp = volumeVariance[0]; 155 for (int j = 1; j <= i; j++) { 156 if (volumeVariance[j] > temp) { 157 temp = volumeVariance[j]; 158 next = j; 159 } 160 } 161 if (temp <= 0.0) { 162 generatedColorCount = i + 1; 163 break; 164 } 165 } 166 167 return new CreateBoxesResult(maxColorCount, generatedColorCount); 168 } 169 createResult(int colorCount)170 List<Integer> createResult(int colorCount) { 171 List<Integer> colors = new ArrayList<>(); 172 for (int i = 0; i < colorCount; ++i) { 173 Box cube = cubes[i]; 174 int weight = volume(cube, weights); 175 if (weight > 0) { 176 int r = volume(cube, momentsR) / weight; 177 int g = volume(cube, momentsG) / weight; 178 int b = volume(cube, momentsB) / weight; 179 int color = (255 << 24) | ((r & 0x0ff) << 16) | ((g & 0x0ff) << 8) | (b & 0x0ff); 180 colors.add(color); 181 } 182 } 183 return colors; 184 } 185 variance(Box cube)186 double variance(Box cube) { 187 int dr = volume(cube, momentsR); 188 int dg = volume(cube, momentsG); 189 int db = volume(cube, momentsB); 190 double xx = 191 moments[getIndex(cube.r1, cube.g1, cube.b1)] 192 - moments[getIndex(cube.r1, cube.g1, cube.b0)] 193 - moments[getIndex(cube.r1, cube.g0, cube.b1)] 194 + moments[getIndex(cube.r1, cube.g0, cube.b0)] 195 - moments[getIndex(cube.r0, cube.g1, cube.b1)] 196 + moments[getIndex(cube.r0, cube.g1, cube.b0)] 197 + moments[getIndex(cube.r0, cube.g0, cube.b1)] 198 - moments[getIndex(cube.r0, cube.g0, cube.b0)]; 199 200 int hypotenuse = dr * dr + dg * dg + db * db; 201 int volume = volume(cube, weights); 202 return xx - hypotenuse / ((double) volume); 203 } 204 cut(Box one, Box two)205 Boolean cut(Box one, Box two) { 206 int wholeR = volume(one, momentsR); 207 int wholeG = volume(one, momentsG); 208 int wholeB = volume(one, momentsB); 209 int wholeW = volume(one, weights); 210 211 MaximizeResult maxRResult = 212 maximize(one, Direction.RED, one.r0 + 1, one.r1, wholeR, wholeG, wholeB, wholeW); 213 MaximizeResult maxGResult = 214 maximize(one, Direction.GREEN, one.g0 + 1, one.g1, wholeR, wholeG, wholeB, wholeW); 215 MaximizeResult maxBResult = 216 maximize(one, Direction.BLUE, one.b0 + 1, one.b1, wholeR, wholeG, wholeB, wholeW); 217 Direction cutDirection; 218 double maxR = maxRResult.maximum; 219 double maxG = maxGResult.maximum; 220 double maxB = maxBResult.maximum; 221 if (maxR >= maxG && maxR >= maxB) { 222 if (maxRResult.cutLocation < 0) { 223 return false; 224 } 225 cutDirection = Direction.RED; 226 } else if (maxG >= maxR && maxG >= maxB) { 227 cutDirection = Direction.GREEN; 228 } else { 229 cutDirection = Direction.BLUE; 230 } 231 232 two.r1 = one.r1; 233 two.g1 = one.g1; 234 two.b1 = one.b1; 235 236 switch (cutDirection) { 237 case RED: 238 one.r1 = maxRResult.cutLocation; 239 two.r0 = one.r1; 240 two.g0 = one.g0; 241 two.b0 = one.b0; 242 break; 243 case GREEN: 244 one.g1 = maxGResult.cutLocation; 245 two.r0 = one.r0; 246 two.g0 = one.g1; 247 two.b0 = one.b0; 248 break; 249 case BLUE: 250 one.b1 = maxBResult.cutLocation; 251 two.r0 = one.r0; 252 two.g0 = one.g0; 253 two.b0 = one.b1; 254 break; 255 } 256 257 one.vol = (one.r1 - one.r0) * (one.g1 - one.g0) * (one.b1 - one.b0); 258 two.vol = (two.r1 - two.r0) * (two.g1 - two.g0) * (two.b1 - two.b0); 259 260 return true; 261 } 262 maximize( Box cube, Direction direction, int first, int last, int wholeR, int wholeG, int wholeB, int wholeW)263 MaximizeResult maximize( 264 Box cube, 265 Direction direction, 266 int first, 267 int last, 268 int wholeR, 269 int wholeG, 270 int wholeB, 271 int wholeW) { 272 int bottomR = bottom(cube, direction, momentsR); 273 int bottomG = bottom(cube, direction, momentsG); 274 int bottomB = bottom(cube, direction, momentsB); 275 int bottomW = bottom(cube, direction, weights); 276 277 double max = 0.0; 278 int cut = -1; 279 280 int halfR = 0; 281 int halfG = 0; 282 int halfB = 0; 283 int halfW = 0; 284 for (int i = first; i < last; i++) { 285 halfR = bottomR + top(cube, direction, i, momentsR); 286 halfG = bottomG + top(cube, direction, i, momentsG); 287 halfB = bottomB + top(cube, direction, i, momentsB); 288 halfW = bottomW + top(cube, direction, i, weights); 289 if (halfW == 0) { 290 continue; 291 } 292 293 double tempNumerator = halfR * halfR + halfG * halfG + halfB * halfB; 294 double tempDenominator = halfW; 295 double temp = tempNumerator / tempDenominator; 296 297 halfR = wholeR - halfR; 298 halfG = wholeG - halfG; 299 halfB = wholeB - halfB; 300 halfW = wholeW - halfW; 301 if (halfW == 0) { 302 continue; 303 } 304 305 tempNumerator = halfR * halfR + halfG * halfG + halfB * halfB; 306 tempDenominator = halfW; 307 temp += (tempNumerator / tempDenominator); 308 309 if (temp > max) { 310 max = temp; 311 cut = i; 312 } 313 } 314 return new MaximizeResult(cut, max); 315 } 316 volume(Box cube, int[] moment)317 static int volume(Box cube, int[] moment) { 318 return (moment[getIndex(cube.r1, cube.g1, cube.b1)] 319 - moment[getIndex(cube.r1, cube.g1, cube.b0)] 320 - moment[getIndex(cube.r1, cube.g0, cube.b1)] 321 + moment[getIndex(cube.r1, cube.g0, cube.b0)] 322 - moment[getIndex(cube.r0, cube.g1, cube.b1)] 323 + moment[getIndex(cube.r0, cube.g1, cube.b0)] 324 + moment[getIndex(cube.r0, cube.g0, cube.b1)] 325 - moment[getIndex(cube.r0, cube.g0, cube.b0)]); 326 } 327 bottom(Box cube, Direction direction, int[] moment)328 static int bottom(Box cube, Direction direction, int[] moment) { 329 switch (direction) { 330 case RED: 331 return -moment[getIndex(cube.r0, cube.g1, cube.b1)] 332 + moment[getIndex(cube.r0, cube.g1, cube.b0)] 333 + moment[getIndex(cube.r0, cube.g0, cube.b1)] 334 - moment[getIndex(cube.r0, cube.g0, cube.b0)]; 335 case GREEN: 336 return -moment[getIndex(cube.r1, cube.g0, cube.b1)] 337 + moment[getIndex(cube.r1, cube.g0, cube.b0)] 338 + moment[getIndex(cube.r0, cube.g0, cube.b1)] 339 - moment[getIndex(cube.r0, cube.g0, cube.b0)]; 340 case BLUE: 341 return -moment[getIndex(cube.r1, cube.g1, cube.b0)] 342 + moment[getIndex(cube.r1, cube.g0, cube.b0)] 343 + moment[getIndex(cube.r0, cube.g1, cube.b0)] 344 - moment[getIndex(cube.r0, cube.g0, cube.b0)]; 345 } 346 throw new IllegalArgumentException("unexpected direction " + direction); 347 } 348 top(Box cube, Direction direction, int position, int[] moment)349 static int top(Box cube, Direction direction, int position, int[] moment) { 350 switch (direction) { 351 case RED: 352 return (moment[getIndex(position, cube.g1, cube.b1)] 353 - moment[getIndex(position, cube.g1, cube.b0)] 354 - moment[getIndex(position, cube.g0, cube.b1)] 355 + moment[getIndex(position, cube.g0, cube.b0)]); 356 case GREEN: 357 return (moment[getIndex(cube.r1, position, cube.b1)] 358 - moment[getIndex(cube.r1, position, cube.b0)] 359 - moment[getIndex(cube.r0, position, cube.b1)] 360 + moment[getIndex(cube.r0, position, cube.b0)]); 361 case BLUE: 362 return (moment[getIndex(cube.r1, cube.g1, position)] 363 - moment[getIndex(cube.r1, cube.g0, position)] 364 - moment[getIndex(cube.r0, cube.g1, position)] 365 + moment[getIndex(cube.r0, cube.g0, position)]); 366 } 367 throw new IllegalArgumentException("unexpected direction " + direction); 368 } 369 370 private static enum Direction { 371 RED, 372 GREEN, 373 BLUE 374 } 375 376 private static final class MaximizeResult { 377 // < 0 if cut impossible 378 int cutLocation; 379 double maximum; 380 MaximizeResult(int cut, double max)381 MaximizeResult(int cut, double max) { 382 this.cutLocation = cut; 383 this.maximum = max; 384 } 385 } 386 387 private static final class CreateBoxesResult { 388 int requestedCount; 389 int resultCount; 390 CreateBoxesResult(int requestedCount, int resultCount)391 CreateBoxesResult(int requestedCount, int resultCount) { 392 this.requestedCount = requestedCount; 393 this.resultCount = resultCount; 394 } 395 } 396 397 private static final class Box { 398 int r0 = 0; 399 int r1 = 0; 400 int g0 = 0; 401 int g1 = 0; 402 int b0 = 0; 403 int b1 = 0; 404 int vol = 0; 405 } 406 } 407