src/trace/filterset.cpp

   1 /*
   2  * Some filters for Potrace in Inkscape
   3  *
   4  * Authors:
   5  *   Bob Jamison <rjamison@titan.com>
   6  *
   7  * Copyright (C) 2004 Bob Jamison
   8  *
   9  * Released under GNU GPL, read the file 'COPYING' for more information
  10  */
  11
  12 #include <stdlib.h>
  13
  14 #include "imagemap-gdk.h"
  15 #include "filterset.h"
  16
  17
  18 /*#########################################################################
  19 ### G A U S S I A N  (smoothing)
  20 #########################################################################*/
  21
  22 /**
  23  *
  24  */
  25 static int gaussMatrix[] =
  26 {
  27      2,  4,  5,  4, 2,
  28      4,  9, 12,  9, 4,
  29      5, 12, 15, 12, 5,
  30      4,  9, 12,  9, 4,
  31      2,  4,  5,  4, 2
  32 };
  33
  34
  35 /**
  36  *
  37  */
  38 GrayMap *grayMapGaussian(GrayMap *me)
  39 {
  40     int width  = me->width;
  41     int height = me->height;
  42     int firstX = 2;
  43     int lastX  = width-3;
  44     int firstY = 2;
  45     int lastY  = height-3;
  46
  47     GrayMap *newGm = GrayMapCreate(width, height);
  48     if (!newGm)
  49         return NULL;
  50
  51     for (int y = 0 ; y<height ; y++)
  52         {
  53         for (int x = 0 ; x<width ; x++)
  54             {
  55             /* image boundaries */
  56             if (x<firstX || x>lastX || y<firstY || y>lastY)
  57                 {
  58                 newGm->setPixel(newGm, x, y, me->getPixel(me, x, y));
  59                 continue;
  60                 }
  61
  62             /* all other pixels */
  63             int gaussIndex = 0;
  64             unsigned long sum = 0;
  65             for (int i= y-2 ; i<=y+2 ; i++)
  66                 {
  67                 for (int j= x-2; j<=x+2 ; j++)
  68                     {
  69                     int weight = gaussMatrix[gaussIndex++];
  70                     sum += me->getPixel(me, j, i) * weight;
  71                     }
  72                 }
  73             sum /= 159;
  74             newGm->setPixel(newGm, x, y, sum);
  75             }
  76         }
  77
  78     return newGm;
  79 }
  80
  81
  82
  83
  84
  85 /**
  86  *
  87  */
  88 RgbMap *rgbMapGaussian(RgbMap *me)
  89 {
  90     int width  = me->width;
  91     int height = me->height;
  92     int firstX = 2;
  93     int lastX  = width-3;
  94     int firstY = 2;
  95     int lastY  = height-3;
  96
  97     RgbMap *newGm = RgbMapCreate(width, height);
  98     if (!newGm)
  99         return NULL;
 100
 101     for (int y = 0 ; y<height ; y++)
 102         {
 103         for (int x = 0 ; x<width ; x++)
 104             {
 105             /* image boundaries */
 106             if (x<firstX || x>lastX || y<firstY || y>lastY)
 107                 {
 108                 newGm->setPixelRGB(newGm, x, y, me->getPixel(me, x, y));
 109                 continue;
 110                 }
 111
 112             /* all other pixels */
 113             int gaussIndex = 0;
 114             int sumR       = 0;
 115             int sumG       = 0;
 116             int sumB       = 0;
 117             for (int i= y-2 ; i<=y+2 ; i++)
 118                 {
 119                 for (int j= x-2; j<=x+2 ; j++)
 120                     {
 121                     int weight = gaussMatrix[gaussIndex++];
 122                     RGB rgb = me->getPixel(me, j, i);
 123                     sumR += weight * (int)rgb.r;
 124                     sumG += weight * (int)rgb.g;
 125                     sumB += weight * (int)rgb.b;
 126                     }
 127                 }
 128             RGB rout;
 129             rout.r = ( sumR / 159 ) & 0xff;
 130             rout.g = ( sumG / 159 ) & 0xff;
 131             rout.b = ( sumB / 159 ) & 0xff;
 132             newGm->setPixelRGB(newGm, x, y, rout);
 133             }
 134         }
 135
 136     return newGm;
 137
 138 }
 139
 140
 141
 142
 143 /*#########################################################################
 144 ### C A N N Y    E D G E    D E T E C T I O N
 145 #########################################################################*/
 146
 147
 148 static int sobelX[] =
 149 {
 150     -1,  0,  1 ,
 151     -2,  0,  2 ,
 152     -1,  0,  1
 153 };
 154
 155 static int sobelY[] =
 156 {
 157      1,  2,  1 ,
 158      0,  0,  0 ,
 159     -1, -2, -1
 160 };
 161
 162
 163
 164 /**
 165  * Perform Sobel convolution on a GrayMap
 166  */
 167 static GrayMap *grayMapSobel(GrayMap *gm,
 168                double dLowThreshold, double dHighThreshold)
 169 {
 170     int width  = gm->width;
 171     int height = gm->height;
 172     int firstX = 1;
 173     int lastX  = width-2;
 174     int firstY = 1;
 175     int lastY  = height-2;
 176
 177     GrayMap *newGm = GrayMapCreate(width, height);
 178     if (!newGm)
 179         return NULL;
 180
 181     for (int y = 0 ; y<height ; y++)
 182         {
 183         for (int x = 0 ; x<width ; x++)
 184             {
 185             unsigned long sum = 0;
 186             /* image boundaries */
 187             if (x<firstX || x>lastX || y<firstY || y>lastY)
 188                 {
 189                 sum = 0;
 190                 }
 191             else
 192                 {
 193                 /* ### SOBEL FILTERING #### */
 194                 long sumX = 0;
 195                 long sumY = 0;
 196                 int sobelIndex = 0;
 197                 for (int i= y-1 ; i<=y+1 ; i++)
 198                     {
 199                     for (int j= x-1; j<=x+1 ; j++)
 200                         {
 201                         sumX += gm->getPixel(gm, j, i) *
 202                              sobelX[sobelIndex++];
 203                         }
 204                     }
 205
 206                 sobelIndex = 0;
 207                 for (int i= y-1 ; i<=y+1 ; i++)
 208                     {
 209                     for (int j= x-1; j<=x+1 ; j++)
 210                         {
 211                         sumY += gm->getPixel(gm, j, i) *
 212                              sobelY[sobelIndex++];
 213                         }
 214                     }
 215                 /*###  GET VALUE ### */
 216                 sum = abs(sumX) + abs(sumY);
 217
 218                 if (sum > 765)
 219                     sum = 765;
 220
 221 #if 0
 222                 /*###  GET ORIENTATION (slow, pedantic way) ### */
 223                 double orient = 0.0;
 224                 if (sumX==0)
 225                     {
 226                     if (sumY==0)
 227                         orient = 0.0;
 228                     else if (sumY<0)
 229                         {
 230                         sumY = -sumY;
 231                         orient = 90.0;
 232                         }
 233                     else
 234                         orient = 90.0;
 235                     }
 236                 else
 237                     {
 238                     orient = 57.295779515 * atan2( ((double)sumY),((double)sumX) );
 239                     if (orient < 0.0)
 240                         orient += 180.0;
 241                     }
 242
 243                 /*###  GET EDGE DIRECTION ### */
 244                 int edgeDirection = 0;
 245                 if (orient < 22.5)
 246                     edgeDirection = 0;
 247                 else if (orient < 67.5)
 248                     edgeDirection = 45;
 249                 else if (orient < 112.5)
 250                     edgeDirection = 90;
 251                 else if (orient < 157.5)
 252                     edgeDirection = 135;
 253 #else
 254                 /*###  GET EDGE DIRECTION (fast way) ### */
 255                 int edgeDirection = 0; /*x,y=0*/
 256                 if (sumX==0)
 257                     {
 258                     if (sumY!=0)
 259                         edgeDirection = 90;
 260                     }
 261                 else
 262                    {
 263                    /*long slope = sumY*1024/sumX;*/
 264                    long slope = (sumY << 10)/sumX;
 265                    if (slope > 2472 || slope< -2472)  /*tan(67.5)*1024*/
 266                        edgeDirection = 90;
 267                    else if (slope > 414) /*tan(22.5)*1024*/
 268                        edgeDirection = 45;
 269                    else if (slope < -414) /*-tan(22.5)*1024*/
 270                        edgeDirection = 135;
 271                    }
 272
 273 #endif
 274                 /* printf("%ld %ld %f %d\n", sumX, sumY, orient, edgeDirection); */
 275
 276                 /*### Get two adjacent pixels in edge direction ### */
 277                 unsigned long leftPixel;
 278                 unsigned long rightPixel;
 279                 if (edgeDirection == 0)
 280                     {
 281                     leftPixel  = gm->getPixel(gm, x-1, y);
 282                     rightPixel = gm->getPixel(gm, x+1, y);
 283                     }
 284                 else if (edgeDirection == 45)
 285                     {
 286                     leftPixel  = gm->getPixel(gm, x-1, y+1);
 287                     rightPixel = gm->getPixel(gm, x+1, y-1);
 288                     }
 289                 else if (edgeDirection == 90)
 290                     {
 291                     leftPixel  = gm->getPixel(gm, x, y-1);
 292                     rightPixel = gm->getPixel(gm, x, y+1);
 293                     }
 294                 else /*135 */
 295                     {
 296                     leftPixel  = gm->getPixel(gm, x-1, y-1);
 297                     rightPixel = gm->getPixel(gm, x+1, y+1);
 298                     }
 299
 300                 /*### Compare current value to adjacent pixels ### */
 301                 /*### if less that either, suppress it ### */
 302                 if (sum < leftPixel || sum < rightPixel)
 303                     sum = 0;
 304                 else
 305                     {
 306                     unsigned long highThreshold =
 307                           (unsigned long)(dHighThreshold * 765.0);
 308                     unsigned long lowThreshold =
 309                           (unsigned long)(dLowThreshold * 765.0);
 310                     if (sum >= highThreshold)
 311                         sum = 765; /* EDGE.  3*255 this needs to be settable */
 312                     else if (sum < lowThreshold)
 313                         sum = 0; /* NONEDGE */
 314                     else
 315                         {
 316                         if ( gm->getPixel(gm, x-1, y-1)> highThreshold ||
 317                              gm->getPixel(gm, x  , y-1)> highThreshold ||
 318                              gm->getPixel(gm, x+1, y-1)> highThreshold ||
 319                              gm->getPixel(gm, x-1, y  )> highThreshold ||
 320                              gm->getPixel(gm, x+1, y  )> highThreshold ||
 321                              gm->getPixel(gm, x-1, y+1)> highThreshold ||
 322                              gm->getPixel(gm, x  , y+1)> highThreshold ||
 323                              gm->getPixel(gm, x+1, y+1)> highThreshold)
 324                             sum = 765; /* EDGE fix me too */
 325                         else
 326                             sum = 0; /* NONEDGE */
 327                         }
 328                     }
 329
 330
 331                 }/* else */
 332             if (sum==0) /* invert light & dark */
 333                 sum = 765;
 334             else
 335                 sum = 0;
 336             newGm->setPixel(newGm, x, y, sum);
 337             }/* for (x) */
 338         }/* for (y) */
 339
 340     return newGm;
 341 }
 342
 343
 344
 345
 346
 347 /**
 348  *
 349  */
 350 GrayMap *
 351 grayMapCanny(GrayMap *gm, double lowThreshold, double highThreshold)
 352 {
 353     if (!gm)
 354         return NULL;
 355     GrayMap *gaussGm = grayMapGaussian(gm);
 356     if (!gaussGm)
 357         return NULL;
 358     /*gaussGm->writePPM(gaussGm, "gauss.ppm");*/
 359
 360     GrayMap *cannyGm = grayMapSobel(gaussGm, lowThreshold, highThreshold);
 361     if (!cannyGm)
 362         return NULL;
 363     /*cannyGm->writePPM(cannyGm, "canny.ppm");*/
 364
 365     gaussGm->destroy(gaussGm);
 366
 367     return cannyGm;
 368 }
 369
 370
 371
 372
 373
 374
 375
 376 /**
 377  *
 378  */
 379 GdkPixbuf *
 380 gdkCanny(GdkPixbuf *img, double lowThreshold, double highThreshold)
 381 {
 382     if (!img)
 383         return NULL;
 384
 385
 386     GrayMap *grayMap = gdkPixbufToGrayMap(img);
 387     if (!grayMap)
 388         return NULL;
 389
 390     /*grayMap->writePPM(grayMap, "gbefore.ppm");*/
 391
 392     GrayMap *cannyGm = grayMapCanny(grayMap,lowThreshold, highThreshold);
 393
 394     grayMap->destroy(grayMap);
 395
 396     if (!cannyGm)
 397         return NULL;
 398
 399     /*grayMap->writePPM(grayMap, "gafter.ppm");*/
 400
 401     GdkPixbuf *newImg = grayMapToGdkPixbuf(cannyGm);
 402
 403
 404     return newImg;
 405 }
 406
 407
 408
 409
 410 /*#########################################################################
 411 ### Q U A N T I Z A T I O N
 412 #########################################################################*/
 413 typedef struct OctreeNode_def OctreeNode;
 414
 415 /**
 416  * The basic octree node
 417  */
 418 struct OctreeNode_def
 419 {
 420     unsigned long r;
 421     unsigned long g;
 422     unsigned long b;
 423     unsigned int  index;
 424     unsigned long nrPixels;
 425     unsigned int  nrChildren;
 426     OctreeNode *parent;
 427     OctreeNode *children[8];
 428 };
 429
 430
 431 /**
 432  * create an octree node, and initialize it
 433  */
 434 OctreeNode *octreeNodeCreate()
 435 {
 436     OctreeNode *node = (OctreeNode *)malloc(sizeof(OctreeNode));
 437     if (!node)
 438         return NULL;
 439     node->r             = 0;
 440     node->g             = 0;
 441     node->b             = 0;
 442     node->index         = 0;
 443     node->nrPixels      = 0;
 444     node->nrChildren    = 0;
 445     node->parent        = NULL;
 446     for (int i=0 ; i<8 ; i++)
 447         node->children[i] = NULL;
 448     return node;
 449 }
 450
 451 /**
 452  *  delete an octree node and its children
 453  */
 454 void octreeNodeDelete(OctreeNode *node)
 455 {
 456     if (!node)
 457         return;
 458     for (int i=0 ; i<8 ; i++)
 459         octreeNodeDelete(node->children[i]);
 460     free(node);
 461 }
 462
 463
 464 /**
 465  *  delete the children of an octree node
 466  */
 467 void octreeNodeDeleteChildren(OctreeNode *node)
 468 {
 469     if (!node)
 470         return;
 471     node->nrChildren = 0;
 472     for (int i=0 ; i<8 ; i++)
 473         {
 474         octreeNodeDelete(node->children[i]);
 475         node->children[i]=NULL;
 476         }
 477 }
 478
 479
 480
 481
 482 /**
 483  *  insert an RGB value into an octree node according to its
 484  *  high-order rgb vector bits
 485  */
 486 int octreeNodeInsert(OctreeNode *root, RGB rgb, int bitsPerSample)
 487 {
 488     OctreeNode *node = root;
 489     int newColor     = FALSE;
 490     int r            = rgb.r;
 491     int g            = rgb.g;
 492     int b            = rgb.b;
 493
 494     int shift = 7;
 495     for (int bit=0 ; bit<bitsPerSample ; bit++)
 496         {
 497         /* update values of all nodes from the root to the leaf */
 498         node->r += r;
 499         node->g += g;
 500         node->b += b;
 501         node->nrPixels++;
 502         int index = (((r >> shift) & 1) << 2) |
 503                     (((g >> shift) & 1) << 1) |
 504                     (((b >> shift) & 1)     ) ;
 505
 506         OctreeNode *child = node->children[index];
 507         if (!child)
 508             {
 509             child                 = octreeNodeCreate();
 510             node->children[index] = child;
 511             child->parent         = node;
 512             node->nrChildren++;
 513             newColor              = TRUE;
 514             }
 515         node = child; /*next level*/
 516         shift--;
 517         }
 518     return newColor;
 519 }
 520
 521
 522
 523
 524
 525 /**
 526  *  find an exact match for an RGB value, at the given bits
 527  *  per sample.  if not found, return -1
 528  */
 529 int octreeNodeFind(OctreeNode *root, RGB rgb, int bitsPerSample)
 530 {
 531     OctreeNode *node = root;
 532     int r            = rgb.r;
 533     int g            = rgb.g;
 534     int b            = rgb.b;
 535
 536     int shift = 7;
 537     for (int bit=0 ; bit<bitsPerSample ; bit++)
 538         {
 539         int index = (((r >> shift) & 1) << 2) |
 540                     (((g >> shift) & 1) << 1) |
 541                     (((b >> shift) & 1)     ) ;
 542
 543         OctreeNode *child = node->children[index];
 544         if (!child)
 545             return -1;
 546         node = child; /*next level*/
 547         shift--;
 548         }
 549     printf("error.  this should not happen\n");
 550     return -1;
 551 }
 552
 553
 554
 555 static void spaces(int nr)
 556 {
 557     for (int i=0; i<nr ; i++)
 558         printf(" ");
 559 }
 560
 561 /**
 562  *  pretty-print an octree node and its children
 563  */
 564 void octreeNodePrint(OctreeNode *node, int indent)
 565 {
 566     spaces(indent); printf("####Node###\n");
 567     spaces(indent); printf("r :%lu\n", node->r);
 568     spaces(indent); printf("g :%lu\n", node->g);
 569     spaces(indent); printf("b :%lu\n", node->b);
 570     spaces(indent); printf("i :%d\n", node->index);
 571     for (unsigned int i=0; i<8; i++)
 572         {
 573         OctreeNode *child = node->children[i];
 574         if (!child)
 575             continue;
 576         spaces(indent); printf("   child %d :", i);
 577         octreeNodePrint(child, indent+4);
 578         }
 579 }
 580
 581 /**
 582  *  Count how many nodes in this (sub)tree
 583  */
 584 static int octreeNodeCount(OctreeNode *node)
 585 {
 586     int count = 1;
 587     for (unsigned int i=0; i<8; i++)
 588         {
 589         OctreeNode *child = node->children[i];
 590         if (!child)
 591             continue;
 592         count += octreeNodeCount(child);
 593         }
 594     return count;
 595 }
 596
 597
 598
 599
 600 /**
 601  * Count all of the leaf nodes in the octree
 602  */
 603 static void octreeLeafArray(OctreeNode *node, OctreeNode **array, int arraySize, int *len)
 604 {
 605     if (!node)
 606         return;
 607     if (node->nrChildren == 0 && *len < arraySize)
 608         {
 609         array[*len] = node;
 610         *len = *len + 1;
 611         }
 612     for (int i=0 ; i<8 ; i++)
 613         octreeLeafArray(node->children[i], array, arraySize, len);
 614 }
 615
 616
 617
 618 /**
 619  *  Count all of the leaf nodes in the octree
 620  */
 621 static int octreeLeafCount(OctreeNode *node)
 622 {
 623     if (!node)
 624         return 0;
 625     if (node->nrChildren == 0)
 626         return 1;
 627     int leaves = 0;
 628     for (int i=0 ; i<8 ; i++)
 629         leaves += octreeLeafCount(node->children[i]);
 630     return leaves;
 631 }
 632
 633 /**
 634  * Mark all of the leaf nodes in the octree with an index nr
 635  */
 636 static void octreeLeafIndex(OctreeNode *node, int *index)
 637 {
 638     if (!node)
 639         return;
 640     if (node->nrChildren == 0)
 641         {
 642         node->index = *index;
 643         *index = *index + 1;
 644         return;
 645         }
 646     for (int i=0 ; i<8 ; i++)
 647         octreeLeafIndex(node->children[i], index);
 648 }
 649
 650
 651
 652 /**
 653  * Find a node that has children, and that also
 654  * has the lowest pixel count
 655  */
 656 static void octreefindLowestLeaf(OctreeNode *node, OctreeNode **lowestLeaf)
 657 {
 658     if (!node)
 659         return;
 660     if (node->nrChildren == 0)
 661         {
 662         if (node->nrPixels < (*lowestLeaf)->nrPixels)
 663             *lowestLeaf = node;
 664         return;
 665         }
 666
 667     for (int i=0 ; i<8 ; i++)
 668         octreefindLowestLeaf(node->children[i], lowestLeaf);
 669 }
 670
 671
 672 /**
 673  * reduce the leaves on an octree to a given number
 674  */
 675 int octreePrune(OctreeNode *root, int nrColors)
 676 {
 677     int leafCount = octreeLeafCount(root);
 678
 679     while (leafCount > nrColors)
 680         {
 681         OctreeNode *lowestLeaf = root;
 682         octreefindLowestLeaf(root, &lowestLeaf);
 683
 684         if (!lowestLeaf)
 685             break; //should never happen
 686
 687        if (lowestLeaf==root)
 688             {
 689             printf("found no leaves\n");
 690             continue;
 691             }
 692
 693         OctreeNode *parent = lowestLeaf->parent;
 694         if (!parent)
 695             continue;
 696
 697         for (int i=0 ; i<8 ; i++)
 698             {
 699             OctreeNode *child = parent->children[i];
 700             if (child == lowestLeaf)
 701                 {
 702                 parent->nrChildren--;
 703                 octreeNodeDelete(child);
 704                 parent->children[i] = NULL;
 705                 break;
 706                 }
 707             }
 708         /*printf("leafCount:%d lowPixels:%lu\n",
 709                leafCount, lowestLeaf->nrPixels);*/
 710         leafCount = octreeLeafCount(root);
 711         }
 712     int index = 0;
 713     octreeLeafIndex(root, &index);
 714
 715     //printf("leafCount:%d\n", leafCount);
 716     //octreeNodePrint(root, 0);
 717
 718     return leafCount;
 719 }
 720
 721
 722
 723 /**
 724  *
 725  */
 726 OctreeNode *octreeBuild(RgbMap *rgbMap, int bitsPerSample, int nrColors)
 727 {
 728     OctreeNode *root = octreeNodeCreate();
 729     if (!root)
 730         return NULL;
 731     for (int y=0 ; y<rgbMap->height ; y++)
 732         {
 733         for (int x=0 ; x<rgbMap->width ; x++)
 734             {
 735             RGB rgb = rgbMap->getPixel(rgbMap, x, y);
 736             octreeNodeInsert(root, rgb, bitsPerSample);
 737             }
 738         }
 739
 740     if (octreePrune(root, nrColors)<0)
 741         {
 742         octreeNodeDelete(root);
 743         return NULL;
 744         }
 745
 746     /* octreeNodePrint(root, 0); */
 747
 748     return root;
 749 }
 750
 751
 752 /* Compare two rgb's for brightness, for the qsort() below */
 753 static int compRGB(const void *a, const void *b)
 754 {
 755     RGB *ra = (RGB *)a;
 756     RGB *rb = (RGB *)b;
 757     int ba = ra->r + ra->g + ra->b;
 758     int bb = rb->r + rb->g + rb->b;
 759     int comp = ba - bb;
 760     return comp;
 761 }
 762
 763 /**
 764  *
 765  */
 766 RGB *makeRGBPalette(OctreeNode *root, int nrColors)
 767 {
 768
 769     OctreeNode **palette = (OctreeNode **)malloc(nrColors * sizeof(OctreeNode *));
 770     if (!palette)
 771         {
 772         return NULL;
 773         }
 774     int len = 0;
 775     octreeLeafArray(root, palette, nrColors, &len);
 776
 777     RGB *rgbpal = (RGB *)malloc(len * sizeof(RGB));
 778     if (!rgbpal)
 779         {
 780         free(palette);
 781         return NULL;
 782         }
 783
 784     for (int i=0; i<len ; i++)
 785         {
 786         OctreeNode *node = palette[i];
 787         RGB rgb;
 788         if (node->nrPixels == 0)
 789             {
 790             continue;
 791             }
 792         rgb.r = (unsigned char)( (node->r / node->nrPixels) & 0xff);
 793         rgb.g = (unsigned char)( (node->g / node->nrPixels) & 0xff);
 794         rgb.b = (unsigned char)( (node->b / node->nrPixels) & 0xff);
 795         int index = node->index;
 796         //printf("Pal: %d %d %d %d\n", rgb.r, rgb.g, rgb.b, index);
 797         rgbpal[index]=rgb;
 798         }
 799
 800     free(palette);
 801
 802     /* sort by brightness, to avoid high contrasts */
 803     qsort((void *)rgbpal, len, sizeof(RGB), compRGB);
 804
 805     return rgbpal;
 806 }
 807
 808
 809 /**
 810  *  Return the closest color in the palette to the request
 811  */
 812 RGB lookupQuantizedRGB(RGB *rgbpalette, int paletteSize, RGB candidate, int *closestIndex)
 813 {
 814     /* slow method */
 815     unsigned long closestMatch = 10000000;
 816     RGB closestRGB = { 0 , 0, 0 };
 817     *closestIndex = 0;
 818     for (int i=0 ; i<paletteSize ; i++)
 819         {
 820         RGB entry = rgbpalette[i];
 821         unsigned long dr    = candidate.r - entry.r;
 822         unsigned long dg    = candidate.g - entry.g;
 823         unsigned long db    = candidate.b - entry.b;
 824         unsigned long match = dr * dr + dg * dg + db * db;
 825         if (match < closestMatch)
 826             {
 827             *closestIndex = i;
 828             closestMatch  = match;
 829             closestRGB    = entry;
 830             }
 831         }
 832
 833     return closestRGB;
 834 }
 835
 836
 837 /**
 838  * Quantize an RGB image to a reduced number of colors.  bitsPerSample
 839  * is usually 3 - 5 out of 8 to conserve cpu and memory
 840  */
 841 IndexedMap *rgbMapQuantize(RgbMap *rgbMap, int bitsPerSample, int nrColors)
 842 {
 843     if (!rgbMap)
 844         return NULL;
 845
 846     OctreeNode *otree = octreeBuild(rgbMap, bitsPerSample, nrColors);
 847     if (!otree)
 848         {
 849         return NULL;
 850         }
 851
 852     /*Make sure we don't request more colors than actually exist*/
 853     int nodeCount = octreeLeafCount(otree);
 854     if (nodeCount < nrColors)
 855         nrColors = nodeCount;
 856
 857     RGB *rgbpal = makeRGBPalette(otree, nrColors);
 858     if (!rgbpal)
 859         {
 860         octreeNodeDelete(otree);
 861         return NULL;
 862         }
 863
 864     /*We have our original and palette. Make the new one*/
 865     IndexedMap *newMap = IndexedMapCreate(rgbMap->width, rgbMap->height);
 866     if (!newMap)
 867         {
 868         free(rgbpal);
 869         octreeNodeDelete(otree);
 870         return NULL;
 871         }
 872
 873     /* fill in the color lookup table */
 874     for (int i=0 ; i< nrColors ; i++)
 875         {
 876         newMap->clut[i] = rgbpal[i];
 877         }
 878     newMap->nrColors = nrColors;
 879
 880     for (int y=0 ; y<rgbMap->height ; y++)
 881         {
 882         for (int x=0 ; x<rgbMap->width ; x++)
 883             {
 884             RGB rgb = rgbMap->getPixel(rgbMap, x, y);
 885             //int indexNr = octreeNodeFind(otree, rgb, bitsPerSample);
 886             //printf("i:%d\n", indexNr);
 887             //RGB quantRgb = rgbpal[indexNr];
 888             int closestIndex;
 889             RGB quantRgb = lookupQuantizedRGB(rgbpal, nrColors, rgb, &closestIndex);
 890             newMap->setPixel(newMap, x, y, (unsigned int)closestIndex);
 891             }
 892         }
 893
 894     free(rgbpal);
 895     octreeNodeDelete(otree);
 896
 897     return newMap;
 898 }
 899
 900
 901
 902 /**
 903  *  Experimental.  Work on this later
 904  */
 905 GrayMap *quantizeBand(RgbMap *rgbMap, int nrColors)
 906 {
 907
 908     int bitsPerSample = 4;
 909
 910     RgbMap *gaussMap = rgbMapGaussian(rgbMap);
 911     //gaussMap->writePPM(gaussMap, "rgbgauss.ppm");
 912
 913     IndexedMap *qMap = rgbMapQuantize(gaussMap, bitsPerSample, nrColors);
 914     //qMap->writePPM(qMap, "rgbquant.ppm");
 915     gaussMap->destroy(gaussMap);
 916
 917     GrayMap *gm = GrayMapCreate(rgbMap->width, rgbMap->height);
 918
 919     /* RGB is quantized.  There should now be a small set of (R+G+B) */
 920     for (int y=0 ; y<qMap->height ; y++)
 921         {
 922         for (int x=0 ; x<qMap->width ; x++)
 923             {
 924             RGB rgb = qMap->getPixelValue(qMap, x, y);
 925             int sum = rgb.r + rgb.g + rgb.b;
 926             if (sum & 1)
 927                 sum = 765;
 928             else
 929                 sum = 0;
 930             /*printf("%d %d %d : %d\n", rgb.r, rgb.g, rgb.b, index);*/
 931             gm->setPixel(gm, x, y, sum);
 932             }
 933         }
 934
 935     qMap->destroy(qMap);
 936
 937     return gm;
 938 }
 939
 940
 941
 942
 943
 944
 945
 946
 947 /*#########################################################################
 948 ### E N D    O F    F I L E
 949 #########################################################################*/
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959