src/trace/quantize.cpp

   1 /*
   2  * Quantization for Inkscape
   3  *
   4  * Authors:
   5  *   Stéphane Gimenez <dev@gim.name>
   6  *
   7  * Copyright (C) 2006 Authors
   8  *
   9  * Released under GNU GPL, read the file 'COPYING' for more information
  10  */
  11
  12 #include <cassert>
  13 #include <cstdio>
  14 #include <stdlib.h>
  15 #include <new>
  16
  17 #include "pool.h"
  18 #include "imagemap.h"
  19
  20 typedef struct Ocnode_def Ocnode;
  21
  22 /**
  23  * an octree node datastructure
  24  */
  25 struct Ocnode_def
  26 {
  27     Ocnode *parent;           // parent node
  28     Ocnode **ref;             // node's reference
  29     Ocnode *child[8];         // children
  30     int nchild;               // number of children
  31     int width;                // width level of this node
  32     RGB rgb;                  // rgb's prefix of that node
  33     unsigned long weight;     // number of pixels this node accounts for
  34     unsigned long rs, gs, bs; // sum of pixels colors this node accounts for
  35     int nleaf;                // number of leaves under this node
  36     unsigned long mi;         // minimum impact
  37 };
  38
  39 /*
  40 -- algorithm principle:
  41
  42 - inspired by the octree method, we associate a tree to a given color map
  43
  44 - nodes in those trees have this shape:
  45
  46                                 parent
  47                                    |
  48         color_prefix(stored in rgb):width
  49      colors_sum(stored in rs,gs,bs)/weight
  50          /               |               \
  51      child1           child2           child3
  52
  53 - (grayscale) trees associated to pixels with colors 87 = 0b1010111 and
  54   69 = 0b1000101 are:
  55
  56            .                 .    <-- roots of the trees
  57            |                 |
  58     1010111:0  and    1000101:0   <-- color prefixes, written in binary form
  59          87/1              69/1   <-- color sums, written in decimal form
  60
  61 - the result of merging the two trees is:
  62
  63                    .
  64                    |
  65                  10:5       <----- longest common prefix and binary width
  66                 156/2       <---.  of the covered color range.
  67             /            \      |
  68     1000101:0      1010111:0    '- sum of colors and quantity of pixels
  69          69/1           87/1       this node accounts for
  70
  71   one should consider three cases when two trees are to be merged:
  72   - one tree range is included in the range of the other one, and the first
  73     tree has to be inserted as a child (or merged with the corresponding
  74     child) of the other.
  75   - their ranges are the same, and their children have to be merged under
  76     a single root.
  77   - ranges have no intersection, and a fork node has to be created (like in
  78     the given example).
  79
  80 - a tree for an image is built dividing the image in 2 parts and merging
  81   the trees obtained recursively for the two parts. a tree for a one pixel
  82   part is a leaf like one of those which were given above.
  83
  84 - last, this tree is reduced a specified number of leaves, deleting first
  85   leaves with minimal impact i.e. [ weight * 2^(2*parentwidth) ] value :
  86   a fair approximation of the impact a leaf removal would have on the final
  87   result : it's the corresponding covered area times the square of the
  88   introduced color distance.
  89
  90   deletion of a node A below a node with only two children is done as
  91   follows :
  92
  93   - when the brother is a leaf, the brother is deleted as well, both nodes
  94     are then represented by their father.
  95
  96      |               |
  97      .       ==>     .
  98     / \
  99    A   .
 100
 101   - otherwise the deletion of A deletes also his father, which plays no
 102     role anymore:
 103
 104      |                |
 105      .       ==>       \
 106     / \                 |
 107    A   .                .
 108       / \              / \
 109
 110   in that way, every leaf removal operation really decreases the remaining
 111   total number of leaves by one.
 112
 113 - very last, color indexes are attributed to leaves; associated colors are
 114   averages, computed from weight and color components sums.
 115
 116 -- improvements to the usual octree method:
 117
 118 - since this algorithm shall often be used to perform quantization using a
 119   very low (2-16) set of colors and not with a usual 256 value, we choose
 120   more carefully which nodes are to be deleted.
 121
 122 - depth of leaves is not fixed to an arbitrary number (which should be 8
 123   when color components are in 0-255), so there is no need to go down to a
 124   depth of 8 for each pixel (at full precision), unless it is really
 125   required.
 126
 127 - tree merging also fastens the overall tree building, and intermediate
 128   processing could be done.
 129
 130 - a huge optimization against the stupid removal algorithm (i.e. find a best
 131   match over the whole tree, remove it and do it again) was implemented:
 132   nodes are marked with the minimal impact of the removal of a leaf below
 133   it. we proceed to the removal recursively. we stop when current removal
 134   level is above the current node minimal, otherwise reached leaves are
 135   removed, and every change over minimal impacts is propagated back to the
 136   whole tree when the recursion ends.
 137
 138 -- specific optimizations
 139
 140 - pool allocation is used to allocate nodes (increased performance on large
 141   images).
 142
 143 */
 144
 145 inline RGB operator>>(RGB rgb, int s)
 146 {
 147   RGB res;
 148   res.r = rgb.r >> s; res.g = rgb.g >> s; res.b = rgb.b >> s;
 149   return res;
 150 }
 151 inline bool operator==(RGB rgb1, RGB rgb2)
 152 {
 153   return (rgb1.r == rgb2.r && rgb1.g == rgb2.g && rgb1.b == rgb2.b);
 154 }
 155
 156 inline int childIndex(RGB rgb)
 157 {
 158     return (((rgb.r)&1)<<2) | (((rgb.g)&1)<<1) | (((rgb.b)&1));
 159 }
 160
 161 /**
 162  * allocate a new node
 163  */
 164 inline Ocnode *ocnodeNew(pool<Ocnode> *pool)
 165 {
 166     Ocnode *node = pool->draw();
 167     node->ref = NULL;
 168     node->parent = NULL;
 169     node->nchild = 0;
 170     for (int i = 0; i < 8; i++) node->child[i] = NULL;
 171     node->mi = 0;
 172     return node;
 173 }
 174
 175 inline void ocnodeFree(pool<Ocnode> *pool, Ocnode *node) {
 176     pool->drop(node);
 177 }
 178
 179
 180 /**
 181  * free a full octree
 182  */
 183 static void octreeDelete(pool<Ocnode> *pool, Ocnode *node)
 184 {
 185     if (!node) return;
 186     for (int i = 0; i < 8; i++)
 187         octreeDelete(pool, node->child[i]);
 188     ocnodeFree(pool, node);
 189 }
 190
 191 /**
 192  *  pretty-print an octree, debugging purposes
 193  */
 194 static void ocnodePrint(Ocnode *node, int indent)
 195 {
 196     if (!node) return;
 197     printf("width:%d weight:%lu rgb:%6x nleaf:%d mi:%lu\n",
 198            node->width,
 199            node->weight,
 200            (unsigned int)(
 201            ((node->rs / node->weight) << 16) +
 202            ((node->gs / node->weight) << 8) +
 203            (node->bs / node->weight)),
 204            node->nleaf,
 205            node->mi
 206            );
 207     for (int i = 0; i < 8; i++) if (node->child[i])
 208         {
 209         for (int k = 0; k < indent; k++) printf(" ");//indentation
 210         printf("[%d:%p] ", i, node->child[i]);
 211         ocnodePrint(node->child[i], indent+2);
 212         }
 213 }
 214 void octreePrint(Ocnode *node)
 215 {
 216     printf("<<octree>>\n");
 217     if (node) printf("[r:%p] ", node); ocnodePrint(node, 2);
 218 }
 219
 220 /**
 221  * builds a single <rgb> color leaf at location <ref>
 222  */
 223 void ocnodeLeaf(pool<Ocnode> *pool, Ocnode **ref, RGB rgb)
 224 {
 225     assert(ref);
 226     Ocnode *node = ocnodeNew(pool);
 227     node->width = 0;
 228     node->rgb = rgb;
 229     node->rs = rgb.r; node->gs = rgb.g; node->bs = rgb.b;
 230     node->weight = 1;
 231     node->nleaf = 1;
 232     node->mi = 0;
 233     node->ref = ref;
 234     *ref = node;
 235 }
 236
 237 /**
 238  *  merge nodes <node1> and <node2> at location <ref> with parent <parent>
 239  */
 240 int octreeMerge(pool<Ocnode> *pool, Ocnode *parent, Ocnode **ref, Ocnode *node1, Ocnode *node2)
 241 {
 242     assert(ref);
 243     if (!node1 && !node2) return 0;
 244     assert(node1 != node2);
 245     if (parent && !*ref) parent->nchild++;
 246     if (!node1)
 247         {
 248         *ref = node2; node2->ref = ref; node2->parent = parent;
 249         return node2->nleaf;
 250         }
 251     if (!node2)
 252         {
 253         *ref = node1; node1->ref = ref; node1->parent = parent;
 254         return node1->nleaf;
 255         }
 256     int dwitdth = node1->width - node2->width;
 257     if (dwitdth > 0 && node1->rgb == node2->rgb >> dwitdth)
 258         {
 259         //place node2 below node1
 260         { *ref = node1; node1->ref = ref; node1->parent = parent; }
 261         int i = childIndex(node2->rgb >> (dwitdth - 1));
 262         node1->rs += node2->rs; node1->gs += node2->gs; node1->bs += node2->bs;
 263         node1->weight += node2->weight;
 264         node1->mi = 0;
 265         if (node1->child[i]) node1->nleaf -= node1->child[i]->nleaf;
 266         node1->nleaf +=
 267           octreeMerge(pool, node1, &node1->child[i], node1->child[i], node2);
 268         return node1->nleaf;
 269         }
 270     else if (dwitdth < 0 && node2->rgb == node1->rgb >> (-dwitdth))
 271         {
 272         //place node1 below node2
 273         { *ref = node2; node2->ref = ref; node2->parent = parent; }
 274         int i = childIndex(node1->rgb >> (-dwitdth - 1));
 275         node2->rs += node1->rs; node2->gs += node1->gs; node2->bs += node1->bs;
 276         node2->weight += node1->weight;
 277         node2->mi = 0;
 278         if (node2->child[i]) node2->nleaf -= node2->child[i]->nleaf;
 279         node2->nleaf +=
 280           octreeMerge(pool, node2, &node2->child[i], node2->child[i], node1);
 281         return node2->nleaf;
 282         }
 283     else
 284         {
 285         //nodes have either no intersection or the same root
 286         Ocnode *newnode;
 287         newnode = ocnodeNew(pool);
 288         newnode->rs = node1->rs + node2->rs;
 289         newnode->gs = node1->gs + node2->gs;
 290         newnode->bs = node1->bs + node2->bs;
 291         newnode->weight = node1->weight + node2->weight;
 292         { *ref = newnode; newnode->ref = ref; newnode->parent = parent; }
 293         if (dwitdth == 0 && node1->rgb == node2->rgb)
 294             {
 295             //merge the nodes in <newnode>
 296             newnode->width = node1->width; // == node2->width
 297             newnode->rgb = node1->rgb;     // == node2->rgb
 298             newnode->nchild = 0;
 299             newnode->nleaf = 0;
 300             if (node1->nchild == 0 && node2->nchild == 0)
 301                 newnode->nleaf = 1;
 302             else
 303                 for (int i = 0; i < 8; i++)
 304                   if (node1->child[i] || node2->child[i])
 305                     newnode->nleaf +=
 306                       octreeMerge(pool, newnode, &newnode->child[i],
 307                                   node1->child[i], node2->child[i]);
 308             ocnodeFree(pool, node1); ocnodeFree(pool, node2);
 309             return newnode->nleaf;
 310             }
 311         else
 312             {
 313             //use <newnode> as a fork node with children <node1> and <node2>
 314             int newwidth =
 315               node1->width > node2->width ? node1->width : node2->width;
 316             RGB rgb1 = node1->rgb >> (newwidth - node1->width);
 317             RGB rgb2 = node2->rgb >> (newwidth - node2->width);
 318             //according to the previous tests <rgb1> != <rgb2> before the loop
 319             while (!(rgb1 == rgb2))
 320               { rgb1 = rgb1 >> 1; rgb2 = rgb2 >> 1; newwidth++; };
 321             newnode->width = newwidth;
 322             newnode->rgb = rgb1;  // == rgb2
 323             newnode->nchild = 2;
 324             newnode->nleaf = node1->nleaf + node2->nleaf;
 325             int i1 = childIndex(node1->rgb >> (newwidth - node1->width - 1));
 326             int i2 = childIndex(node2->rgb >> (newwidth - node2->width - 1));
 327             node1->parent = newnode;
 328             node1->ref = &newnode->child[i1];
 329             newnode->child[i1] = node1;
 330             node2->parent = newnode;
 331             node2->ref = &newnode->child[i2];
 332             newnode->child[i2] = node2;
 333             return newnode->nleaf;
 334             }
 335         }
 336 }
 337
 338 /**
 339  * upatade mi value for leaves
 340  */
 341 static inline void ocnodeMi(Ocnode *node)
 342 {
 343     node->mi = node->parent ?
 344        node->weight << (2 * node->parent->width) : 0;
 345 }
 346
 347 /**
 348  * remove leaves whose prune impact value is lower than <lvl>. at most
 349  * <count> leaves are removed, and <count> is decreased on each removal.
 350  * all parameters including minimal impact values are regenerated.
 351  */
 352 static void ocnodeStrip(pool<Ocnode> *pool, Ocnode **ref, int *count, unsigned long lvl)
 353 {
 354     Ocnode *node = *ref;
 355     if (!count || !node) return;
 356     assert(ref == node->ref);
 357     if (node->nchild == 0) // leaf node
 358         {
 359         if (!node->mi) ocnodeMi(node); //mi generation may be required
 360         if (node->mi > lvl) return; //leaf is above strip level
 361         ocnodeFree(pool, node);
 362         *ref = NULL;
 363         (*count)--;
 364         }
 365     else
 366         {
 367         if (node->mi && node->mi > lvl) //node is above strip level
 368             return;
 369         node->nchild = 0;
 370         node->nleaf = 0;
 371         node->mi = 0;
 372         Ocnode **lonelychild = NULL;
 373         for (int i = 0; i < 8; i++) if (node->child[i])
 374             {
 375             ocnodeStrip(pool, &node->child[i], count, lvl);
 376             if (node->child[i])
 377                 {
 378                 lonelychild = &node->child[i];
 379                 node->nchild++;
 380                 node->nleaf += node->child[i]->nleaf;
 381                 if (!node->mi || node->mi > node->child[i]->mi)
 382                     node->mi = node->child[i]->mi;
 383                 }
 384             }
 385       // tree adjustments
 386         if (node->nchild == 0)
 387             {
 388             (*count)++;
 389             node->nleaf = 1;
 390             ocnodeMi(node);
 391             }
 392         else if (node->nchild == 1)
 393             {
 394             if ((*lonelychild)->nchild == 0)
 395                 {
 396                 //remove the <lonelychild> leaf under a 1 child node
 397                 node->nchild = 0;
 398                 node->nleaf = 1;
 399                 ocnodeMi(node);
 400                 ocnodeFree(pool, *lonelychild);
 401                 *lonelychild = NULL;
 402                 }
 403             else
 404                 {
 405                 //make a bridge to <lonelychild> over a 1 child node
 406                 (*lonelychild)->parent = node->parent;
 407                 (*lonelychild)->ref = ref;
 408                 ocnodeFree(pool, node);
 409                 *ref = *lonelychild;
 410                 }
 411             }
 412         }
 413 }
 414
 415 /**
 416  * reduce the leaves of an octree to a given number
 417  */
 418 void octreePrune(pool<Ocnode> *pool, Ocnode **ref, int ncolor)
 419 {
 420   assert(ref);
 421   assert(ncolor > 0);
 422   //printf("pruning down to %d colors:\n", ncolor);//debug
 423   int n = (*ref)->nleaf - ncolor;
 424   if (!*ref || n <= 0) return;
 425   while (n > 0)
 426       {
 427       //printf("removals to go: %10d\t", n);//debug
 428       //printf("current prune impact: %10lu\n", (*ref)->mi);//debug
 429       //calling strip with global minimum impact of the tree
 430       ocnodeStrip(pool, ref, &n, (*ref)->mi);
 431       }
 432 }
 433
 434 /**
 435  * build an octree associated to the area of a color map <rgbmap>,
 436  * included in the specified (x1,y1)--(x2,y2) rectangle.
 437  */
 438 void octreeBuildArea(pool<Ocnode> *pool, RgbMap *rgbmap, Ocnode **ref,
 439                      int x1, int y1, int x2, int y2, int ncolor)
 440 {
 441     int dx = x2 - x1, dy = y2 - y1;
 442     int xm = x1 + dx/2, ym = y1 + dy/2;
 443     Ocnode *ref1 = NULL;
 444     Ocnode *ref2 = NULL;
 445     if (dx == 1 && dy == 1)
 446         ocnodeLeaf(pool, ref, rgbmap->getPixel(rgbmap, x1, y1));
 447     else if (dx > dy)
 448         {
 449         octreeBuildArea(pool, rgbmap, &ref1, x1, y1, xm, y2, ncolor);
 450         octreeBuildArea(pool, rgbmap, &ref2, xm, y1, x2, y2, ncolor);
 451         octreeMerge(pool, NULL, ref, ref1, ref2);
 452         }
 453     else
 454         {
 455         octreeBuildArea(pool, rgbmap, &ref1, x1, y1, x2, ym, ncolor);
 456         octreeBuildArea(pool, rgbmap, &ref2, x1, ym, x2, y2, ncolor);
 457         octreeMerge(pool, NULL, ref, ref1, ref2);
 458         }
 459
 460     //octreePrune(ref, 2*ncolor);
 461     //affects result quality for almost same performance :/
 462 }
 463
 464 /**
 465  * build an octree associated to the <rgbmap> color map,
 466  * pruned to <ncolor> colors.
 467  */
 468 Ocnode *octreeBuild(pool<Ocnode> *pool, RgbMap *rgbmap, int ncolor)
 469 {
 470     //create the octree
 471     Ocnode *node = NULL;
 472     octreeBuildArea(pool,
 473                     rgbmap, &node,
 474                     0, 0, rgbmap->width, rgbmap->height, ncolor
 475                     );
 476
 477     //prune the octree
 478     octreePrune(pool, &node, ncolor);
 479
 480     //octreePrint(node);//debug
 481
 482     return node;
 483 }
 484
 485 /**
 486  * compute the color palette associated to an octree.
 487  */
 488 static void octreeIndex(Ocnode *node, RGB *rgbpal, int *index)
 489 {
 490     if (!node) return;
 491     if (node->nchild == 0)
 492         {
 493         rgbpal[*index].r = node->rs / node->weight;
 494         rgbpal[*index].g = node->gs / node->weight;
 495         rgbpal[*index].b = node->bs / node->weight;
 496         (*index)++;
 497         }
 498     else
 499         for (int i = 0; i < 8; i++)
 500             if (node->child[i])
 501                 octreeIndex(node->child[i], rgbpal, index);
 502 }
 503
 504 /**
 505  * compute the squared distance between two colors
 506  */
 507 static int distRGB(RGB rgb1, RGB rgb2)
 508 {
 509     return
 510       (rgb1.r - rgb2.r) * (rgb1.r - rgb2.r)
 511     + (rgb1.g - rgb2.g) * (rgb1.g - rgb2.g)
 512     + (rgb1.b - rgb2.b) * (rgb1.b - rgb2.b);
 513 }
 514
 515 /**
 516  * find the index of closest color in a palette
 517  */
 518 static int findRGB(RGB *rgbpal, int ncolor, RGB rgb)
 519 {
 520     //assert(ncolor > 0);
 521     //assert(rgbpal);
 522     int index = -1, dist = 0;
 523     for (int k = 0; k < ncolor; k++)
 524         {
 525         int d = distRGB(rgbpal[k], rgb);
 526         if (index == -1 || d < dist) { dist = d; index = k; }
 527         }
 528     return index;
 529 }
 530
 531 /**
 532  * (qsort) compare two colors for brightness
 533  */
 534 static int compRGB(const void *a, const void *b)
 535 {
 536     RGB *ra = (RGB *)a, *rb = (RGB *)b;
 537     return (ra->r + ra->g + ra->b) - (rb->r + rb->g + rb->b);
 538 }
 539
 540 /**
 541  * quantize an RGB image to a reduced number of colors.
 542  */
 543 IndexedMap *rgbMapQuantize(RgbMap *rgbmap, int ncolor)
 544 {
 545     assert(rgbmap);
 546     assert(ncolor > 0);
 547
 548     pool<Ocnode> pool;
 549
 550     Ocnode *tree;
 551     try {
 552       tree = octreeBuild(&pool, rgbmap, ncolor);
 553     }
 554     catch (std::bad_alloc& ex) {
 555       return NULL; //should do smthg else?
 556     }
 557
 558     RGB *rgbpal = new RGB[ncolor];
 559     int indexes = 0;
 560     octreeIndex(tree, rgbpal, &indexes);
 561
 562     octreeDelete(&pool, tree);
 563
 564     // stacking with increasing contrasts
 565     qsort((void *)rgbpal, indexes, sizeof(RGB), compRGB);
 566
 567     // make the new map
 568     IndexedMap *newmap = IndexedMapCreate(rgbmap->width, rgbmap->height);
 569     if (!newmap) { delete rgbpal; return NULL; }
 570
 571     // fill in the color lookup table
 572     for (int i = 0; i < indexes; i++) newmap->clut[i] = rgbpal[i];
 573     newmap->nrColors = indexes;
 574
 575     // fill in new map pixels
 576     for (int y = 0; y < rgbmap->height; y++)
 577         {
 578         for (int x = 0; x < rgbmap->width; x++)
 579             {
 580             RGB rgb = rgbmap->getPixel(rgbmap, x, y);
 581             int index = findRGB(rgbpal, ncolor, rgb);
 582             newmap->setPixel(newmap, x, y, index);
 583             }
 584         }
 585
 586     delete rgbpal;
 587     return newmap;
 588 }