QuantOctree.c 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
  1. /* Copyright (c) 2010 Oliver Tonnhofer <olt@bogosoft.com>, Omniscale
  2. //
  3. // Permission is hereby granted, free of charge, to any person obtaining a copy
  4. // of this software and associated documentation files (the "Software"), to deal
  5. // in the Software without restriction, including without limitation the rights
  6. // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  7. // copies of the Software, and to permit persons to whom the Software is
  8. // furnished to do so, subject to the following conditions:
  9. //
  10. // The above copyright notice and this permission notice shall be included in
  11. // all copies or substantial portions of the Software.
  12. //
  13. // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  14. // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  15. // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  16. // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  17. // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  18. // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  19. // THE SOFTWARE.
  20. */
  21. /*
  22. // This file implements a variation of the octree color quantization algorithm.
  23. */
  24. #include <stdio.h>
  25. #include <stdlib.h>
  26. #include <string.h>
  27. #include <limits.h>
  28. #include "QuantOctree.h"
  29. typedef struct _ColorBucket{
  30. /* contains palette index when used for look up cube */
  31. uint32_t count;
  32. uint64_t r;
  33. uint64_t g;
  34. uint64_t b;
  35. uint64_t a;
  36. } *ColorBucket;
  37. typedef struct _ColorCube{
  38. unsigned int rBits, gBits, bBits, aBits;
  39. unsigned int rWidth, gWidth, bWidth, aWidth;
  40. unsigned int rOffset, gOffset, bOffset, aOffset;
  41. long size;
  42. ColorBucket buckets;
  43. } *ColorCube;
  44. #define MAX(a, b) (a)>(b) ? (a) : (b)
  45. static ColorCube
  46. new_color_cube(int r, int g, int b, int a) {
  47. ColorCube cube;
  48. /* malloc check ok, small constant allocation */
  49. cube = malloc(sizeof(struct _ColorCube));
  50. if (!cube) return NULL;
  51. cube->rBits = MAX(r, 0);
  52. cube->gBits = MAX(g, 0);
  53. cube->bBits = MAX(b, 0);
  54. cube->aBits = MAX(a, 0);
  55. /* overflow check for size multiplication below */
  56. if (cube->rBits + cube->gBits + cube->bBits + cube->aBits > 31) {
  57. free(cube);
  58. return NULL;
  59. }
  60. /* the width of the cube for each dimension */
  61. cube->rWidth = 1<<cube->rBits;
  62. cube->gWidth = 1<<cube->gBits;
  63. cube->bWidth = 1<<cube->bBits;
  64. cube->aWidth = 1<<cube->aBits;
  65. /* the offsets of each color */
  66. cube->rOffset = cube->gBits + cube->bBits + cube->aBits;
  67. cube->gOffset = cube->bBits + cube->aBits;
  68. cube->bOffset = cube->aBits;
  69. cube->aOffset = 0;
  70. /* the number of color buckets */
  71. cube->size = cube->rWidth * cube->gWidth * cube->bWidth * cube->aWidth;
  72. /* malloc check ok, overflow checked above */
  73. cube->buckets = calloc(cube->size, sizeof(struct _ColorBucket));
  74. if (!cube->buckets) {
  75. free(cube);
  76. return NULL;
  77. }
  78. return cube;
  79. }
  80. static void
  81. free_color_cube(ColorCube cube) {
  82. if (cube != NULL) {
  83. free(cube->buckets);
  84. free(cube);
  85. }
  86. }
  87. static long
  88. color_bucket_offset_pos(const ColorCube cube,
  89. unsigned int r, unsigned int g, unsigned int b, unsigned int a)
  90. {
  91. return r<<cube->rOffset | g<<cube->gOffset | b<<cube->bOffset | a<<cube->aOffset;
  92. }
  93. static long
  94. color_bucket_offset(const ColorCube cube, const Pixel *p) {
  95. unsigned int r = p->c.r>>(8-cube->rBits);
  96. unsigned int g = p->c.g>>(8-cube->gBits);
  97. unsigned int b = p->c.b>>(8-cube->bBits);
  98. unsigned int a = p->c.a>>(8-cube->aBits);
  99. return color_bucket_offset_pos(cube, r, g, b, a);
  100. }
  101. static ColorBucket
  102. color_bucket_from_cube(const ColorCube cube, const Pixel *p) {
  103. unsigned int offset = color_bucket_offset(cube, p);
  104. return &cube->buckets[offset];
  105. }
  106. static void
  107. add_color_to_color_cube(const ColorCube cube, const Pixel *p) {
  108. ColorBucket bucket = color_bucket_from_cube(cube, p);
  109. bucket->count += 1;
  110. bucket->r += p->c.r;
  111. bucket->g += p->c.g;
  112. bucket->b += p->c.b;
  113. bucket->a += p->c.a;
  114. }
  115. static long
  116. count_used_color_buckets(const ColorCube cube) {
  117. long usedBuckets = 0;
  118. long i;
  119. for (i=0; i < cube->size; i++) {
  120. if (cube->buckets[i].count > 0) {
  121. usedBuckets += 1;
  122. }
  123. }
  124. return usedBuckets;
  125. }
  126. static void
  127. avg_color_from_color_bucket(const ColorBucket bucket, Pixel *dst) {
  128. float count = bucket->count;
  129. if (count != 0) {
  130. dst->c.r = (int)(bucket->r / count);
  131. dst->c.g = (int)(bucket->g / count);
  132. dst->c.b = (int)(bucket->b / count);
  133. dst->c.a = (int)(bucket->a / count);
  134. } else {
  135. dst->c.r = 0;
  136. dst->c.g = 0;
  137. dst->c.b = 0;
  138. dst->c.a = 0;
  139. }
  140. }
  141. static int
  142. compare_bucket_count(const ColorBucket a, const ColorBucket b) {
  143. return b->count - a->count;
  144. }
  145. static ColorBucket
  146. create_sorted_color_palette(const ColorCube cube) {
  147. ColorBucket buckets;
  148. if (cube->size > LONG_MAX / sizeof(struct _ColorBucket)) {
  149. return NULL;
  150. }
  151. /* malloc check ok, calloc + overflow check above for memcpy */
  152. buckets = calloc(cube->size, sizeof(struct _ColorBucket));
  153. if (!buckets) return NULL;
  154. memcpy(buckets, cube->buckets, sizeof(struct _ColorBucket)*cube->size);
  155. qsort(buckets, cube->size, sizeof(struct _ColorBucket),
  156. (int (*)(void const *, void const *))&compare_bucket_count);
  157. return buckets;
  158. }
  159. void add_bucket_values(ColorBucket src, ColorBucket dst) {
  160. dst->count += src->count;
  161. dst->r += src->r;
  162. dst->g += src->g;
  163. dst->b += src->b;
  164. dst->a += src->a;
  165. }
  166. /* expand or shrink a given cube to level */
  167. static ColorCube copy_color_cube(const ColorCube cube,
  168. int rBits, int gBits, int bBits, int aBits)
  169. {
  170. unsigned int r, g, b, a;
  171. long src_pos, dst_pos;
  172. unsigned int src_reduce[4] = {0}, dst_reduce[4] = {0};
  173. unsigned int width[4];
  174. ColorCube result;
  175. result = new_color_cube(rBits, gBits, bBits, aBits);
  176. if (!result) return NULL;
  177. if (cube->rBits > rBits) {
  178. dst_reduce[0] = cube->rBits - result->rBits;
  179. width[0] = cube->rWidth;
  180. } else {
  181. src_reduce[0] = result->rBits - cube->rBits;
  182. width[0] = result->rWidth;
  183. }
  184. if (cube->gBits > gBits) {
  185. dst_reduce[1] = cube->gBits - result->gBits;
  186. width[1] = cube->gWidth;
  187. } else {
  188. src_reduce[1] = result->gBits - cube->gBits;
  189. width[1] = result->gWidth;
  190. }
  191. if (cube->bBits > bBits) {
  192. dst_reduce[2] = cube->bBits - result->bBits;
  193. width[2] = cube->bWidth;
  194. } else {
  195. src_reduce[2] = result->bBits - cube->bBits;
  196. width[2] = result->bWidth;
  197. }
  198. if (cube->aBits > aBits) {
  199. dst_reduce[3] = cube->aBits - result->aBits;
  200. width[3] = cube->aWidth;
  201. } else {
  202. src_reduce[3] = result->aBits - cube->aBits;
  203. width[3] = result->aWidth;
  204. }
  205. for (r=0; r<width[0]; r++) {
  206. for (g=0; g<width[1]; g++) {
  207. for (b=0; b<width[2]; b++) {
  208. for (a=0; a<width[3]; a++) {
  209. src_pos = color_bucket_offset_pos(cube,
  210. r>>src_reduce[0],
  211. g>>src_reduce[1],
  212. b>>src_reduce[2],
  213. a>>src_reduce[3]);
  214. dst_pos = color_bucket_offset_pos(result,
  215. r>>dst_reduce[0],
  216. g>>dst_reduce[1],
  217. b>>dst_reduce[2],
  218. a>>dst_reduce[3]);
  219. add_bucket_values(
  220. &cube->buckets[src_pos],
  221. &result->buckets[dst_pos]
  222. );
  223. }
  224. }
  225. }
  226. }
  227. return result;
  228. }
  229. void
  230. subtract_color_buckets(ColorCube cube, ColorBucket buckets, long nBuckets) {
  231. ColorBucket minuend, subtrahend;
  232. long i;
  233. Pixel p;
  234. for (i=0; i<nBuckets; i++) {
  235. subtrahend = &buckets[i];
  236. // If the subtrahend contains no buckets, there is nothing to subtract.
  237. if (subtrahend->count == 0) continue;
  238. avg_color_from_color_bucket(subtrahend, &p);
  239. minuend = color_bucket_from_cube(cube, &p);
  240. minuend->count -= subtrahend->count;
  241. minuend->r -= subtrahend->r;
  242. minuend->g -= subtrahend->g;
  243. minuend->b -= subtrahend->b;
  244. minuend->a -= subtrahend->a;
  245. }
  246. }
  247. static void
  248. set_lookup_value(const ColorCube cube, const Pixel *p, long value) {
  249. ColorBucket bucket = color_bucket_from_cube(cube, p);
  250. bucket->count = value;
  251. }
  252. uint64_t
  253. lookup_color(const ColorCube cube, const Pixel *p) {
  254. ColorBucket bucket = color_bucket_from_cube(cube, p);
  255. return bucket->count;
  256. }
  257. void add_lookup_buckets(ColorCube cube, ColorBucket palette, long nColors, long offset) {
  258. long i;
  259. Pixel p;
  260. for (i=offset; i<offset+nColors; i++) {
  261. avg_color_from_color_bucket(&palette[i], &p);
  262. set_lookup_value(cube, &p, i);
  263. }
  264. }
  265. ColorBucket
  266. combined_palette(ColorBucket bucketsA, long nBucketsA, ColorBucket bucketsB, long nBucketsB) {
  267. ColorBucket result;
  268. if (nBucketsA > LONG_MAX - nBucketsB ||
  269. (nBucketsA+nBucketsB) > LONG_MAX / sizeof(struct _ColorBucket)) {
  270. return NULL;
  271. }
  272. /* malloc check ok, overflow check above */
  273. result = calloc(nBucketsA + nBucketsB, sizeof(struct _ColorBucket));
  274. if (!result) {
  275. return NULL;
  276. }
  277. memcpy(result, bucketsA, sizeof(struct _ColorBucket) * nBucketsA);
  278. memcpy(&result[nBucketsA], bucketsB, sizeof(struct _ColorBucket) * nBucketsB);
  279. return result;
  280. }
  281. static Pixel *
  282. create_palette_array(const ColorBucket palette, unsigned int paletteLength) {
  283. Pixel *paletteArray;
  284. unsigned int i;
  285. /* malloc check ok, calloc for overflow */
  286. paletteArray = calloc(paletteLength, sizeof(Pixel));
  287. if (!paletteArray) return NULL;
  288. for (i=0; i<paletteLength; i++) {
  289. avg_color_from_color_bucket(&palette[i], &paletteArray[i]);
  290. }
  291. return paletteArray;
  292. }
  293. static void
  294. map_image_pixels(const Pixel *pixelData,
  295. uint32_t nPixels,
  296. const ColorCube lookupCube,
  297. uint32_t *pixelArray)
  298. {
  299. long i;
  300. for (i=0; i<nPixels; i++) {
  301. pixelArray[i] = lookup_color(lookupCube, &pixelData[i]);
  302. }
  303. }
  304. const int CUBE_LEVELS[8] = {4, 4, 4, 0, 2, 2, 2, 0};
  305. const int CUBE_LEVELS_ALPHA[8] = {3, 4, 3, 3, 2, 2, 2, 2};
  306. int quantize_octree(Pixel *pixelData,
  307. uint32_t nPixels,
  308. uint32_t nQuantPixels,
  309. Pixel **palette,
  310. uint32_t *paletteLength,
  311. uint32_t **quantizedPixels,
  312. int withAlpha)
  313. {
  314. ColorCube fineCube = NULL;
  315. ColorCube coarseCube = NULL;
  316. ColorCube lookupCube = NULL;
  317. ColorCube coarseLookupCube = NULL;
  318. ColorBucket paletteBucketsCoarse = NULL;
  319. ColorBucket paletteBucketsFine = NULL;
  320. ColorBucket paletteBuckets = NULL;
  321. uint32_t *qp = NULL;
  322. long i;
  323. long nCoarseColors, nFineColors, nAlreadySubtracted;
  324. const int *cubeBits;
  325. if (withAlpha) {
  326. cubeBits = CUBE_LEVELS_ALPHA;
  327. }
  328. else {
  329. cubeBits = CUBE_LEVELS;
  330. }
  331. /*
  332. Create two color cubes, one fine grained with 8x16x8=1024
  333. colors buckets and a coarse with 4x4x4=64 color buckets.
  334. The coarse one guarantees that there are color buckets available for
  335. the whole color range (assuming nQuantPixels > 64).
  336. For a quantization to 256 colors all 64 coarse colors will be used
  337. plus the 192 most used color buckets from the fine color cube.
  338. The average of all colors within one bucket is used as the actual
  339. color for that bucket.
  340. For images with alpha the cubes gets a forth dimension,
  341. 8x16x8x8 and 4x4x4x4.
  342. */
  343. /* create fine cube */
  344. fineCube = new_color_cube(cubeBits[0], cubeBits[1],
  345. cubeBits[2], cubeBits[3]);
  346. if (!fineCube) goto error;
  347. for (i=0; i<nPixels; i++) {
  348. add_color_to_color_cube(fineCube, &pixelData[i]);
  349. }
  350. /* create coarse cube */
  351. coarseCube = copy_color_cube(fineCube, cubeBits[4], cubeBits[5],
  352. cubeBits[6], cubeBits[7]);
  353. if (!coarseCube) goto error;
  354. nCoarseColors = count_used_color_buckets(coarseCube);
  355. /* limit to nQuantPixels */
  356. if (nCoarseColors > nQuantPixels)
  357. nCoarseColors = nQuantPixels;
  358. /* how many space do we have in our palette for fine colors? */
  359. nFineColors = nQuantPixels - nCoarseColors;
  360. /* create fine color palette */
  361. paletteBucketsFine = create_sorted_color_palette(fineCube);
  362. if (!paletteBucketsFine) goto error;
  363. /* remove the used fine colors from the coarse cube */
  364. subtract_color_buckets(coarseCube, paletteBucketsFine, nFineColors);
  365. /* did the subtraction cleared one or more coarse bucket? */
  366. while (nCoarseColors > count_used_color_buckets(coarseCube)) {
  367. /* then we can use the free buckets for fine colors */
  368. nAlreadySubtracted = nFineColors;
  369. nCoarseColors = count_used_color_buckets(coarseCube);
  370. nFineColors = nQuantPixels - nCoarseColors;
  371. subtract_color_buckets(coarseCube, &paletteBucketsFine[nAlreadySubtracted],
  372. nFineColors-nAlreadySubtracted);
  373. }
  374. /* create our palette buckets with fine and coarse combined */
  375. paletteBucketsCoarse = create_sorted_color_palette(coarseCube);
  376. if (!paletteBucketsCoarse) goto error;
  377. paletteBuckets = combined_palette(paletteBucketsCoarse, nCoarseColors,
  378. paletteBucketsFine, nFineColors);
  379. free(paletteBucketsFine);
  380. paletteBucketsFine = NULL;
  381. free(paletteBucketsCoarse);
  382. paletteBucketsCoarse = NULL;
  383. if (!paletteBuckets) goto error;
  384. /* add all coarse colors to our coarse lookup cube. */
  385. coarseLookupCube = new_color_cube(cubeBits[4], cubeBits[5],
  386. cubeBits[6], cubeBits[7]);
  387. if (!coarseLookupCube) goto error;
  388. add_lookup_buckets(coarseLookupCube, paletteBuckets, nCoarseColors, 0);
  389. /* expand coarse cube (64) to larger fine cube (4k). the value of each
  390. coarse bucket is then present in the according 64 fine buckets. */
  391. lookupCube = copy_color_cube(coarseLookupCube, cubeBits[0], cubeBits[1],
  392. cubeBits[2], cubeBits[3]);
  393. if (!lookupCube) goto error;
  394. /* add fine colors to the lookup cube */
  395. add_lookup_buckets(lookupCube, paletteBuckets, nFineColors, nCoarseColors);
  396. /* create result pixels and map palette indices */
  397. /* malloc check ok, calloc for overflow */
  398. qp = calloc(nPixels, sizeof(Pixel));
  399. if (!qp) goto error;
  400. map_image_pixels(pixelData, nPixels, lookupCube, qp);
  401. /* convert palette buckets to RGB pixel palette */
  402. *palette = create_palette_array(paletteBuckets, nQuantPixels);
  403. if (!(*palette)) goto error;
  404. *quantizedPixels = qp;
  405. *paletteLength = nQuantPixels;
  406. free_color_cube(coarseCube);
  407. free_color_cube(fineCube);
  408. free_color_cube(lookupCube);
  409. free_color_cube(coarseLookupCube);
  410. free(paletteBuckets);
  411. return 1;
  412. error:
  413. /* everything is initialized to NULL
  414. so we are safe to call free */
  415. free(qp);
  416. free_color_cube(lookupCube);
  417. free_color_cube(coarseLookupCube);
  418. free(paletteBuckets);
  419. free(paletteBucketsCoarse);
  420. free(paletteBucketsFine);
  421. free_color_cube(coarseCube);
  422. free_color_cube(fineCube);
  423. return 0;
  424. }