src/trace/siox.cpp

   1 /*
   2    Copyright 2005, 2006 by Gerald Friedland, Kristian Jantz and Lars Knipping
   3
   4    Conversion to C++ for Inkscape by Bob Jamison
   5
   6    Licensed under the Apache License, Version 2.0 (the "License");
   7    you may not use this file except in compliance with the License.
   8    You may obtain a copy of the License at
   9
  10    http://www.apache.org/licenses/LICENSE-2.0
  11
  12    Unless required by applicable law or agreed to in writing, software
  13    distributed under the License is distributed on an "AS IS" BASIS,
  14    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  15    See the License for the specific language governing permissions and
  16    limitations under the License.
  17  */
  18 #include "siox.h"
  19
  20 #include <math.h>
  21 #include <stdarg.h>
  22 #include <map>
  23
  24
  25 namespace org
  26 {
  27
  28 namespace siox
  29 {
  30
  31
  32
  33 //########################################################################
  34 //#  C L A B
  35 //########################################################################
  36
  37 static std::map<unsigned long, CLAB> clabLookupTable;
  38
  39 /**
  40  * Convert integer A, R, G, B values into an pixel value.
  41  */
  42 static unsigned long getRGB(int a, int r, int g, int b)
  43 {
  44     if (a<0)  a=0;
  45     else if (a>255) a=255;
  46
  47     if (r<0) r=0;
  48     else if (r>255) r=255;
  49
  50     if (g<0) g=0;
  51     else if (g>255) g=255;
  52
  53     if (b<0) b=0;
  54     else if (b>255) b=255;
  55
  56     return (a<<24)|(r<<16)|(g<<8)|b;
  57 }
  58
  59
  60
  61 /**
  62  * Convert float A, R, G, B values (0.0-1.0) into an pixel value.
  63  */
  64 static unsigned long getRGB(float a, float r, float g, float b)
  65 {
  66     return getRGB((int)(a * 256.0),
  67                   (int)(r * 256.0),
  68                   (int)(g * 256.0),
  69                   (int)(b * 256.0));
  70 }
  71
  72
  73
  74 /**
  75  * Construct this CLAB from a packed-pixel ARGB value
  76  */
  77 CLAB::CLAB(unsigned long rgb)
  78 {
  79     //First try looking up in the cache
  80     std::map<unsigned long, CLAB>::iterator iter;
  81     iter = clabLookupTable.find(rgb);
  82     if (iter != clabLookupTable.end())
  83         {
  84         CLAB res = iter->second;
  85         C = res.C;
  86         L = res.L;
  87         A = res.A;
  88         B = res.B;
  89         }
  90
  91
  92     int ir  = (rgb>>16) & 0xff;
  93     int ig  = (rgb>> 8) & 0xff;
  94     int ib  = (rgb    ) & 0xff;
  95
  96     float fr = ((float)ir) / 255.0;
  97     float fg = ((float)ig) / 255.0;
  98     float fb = ((float)ib) / 255.0;
  99
 100     if (fr > 0.04045)
 101         fr = (float) pow((fr + 0.055) / 1.055, 2.4);
 102     else
 103         fr = fr / 12.92;
 104
 105     if (fg > 0.04045)
 106         fg = (float) pow((fg + 0.055) / 1.055, 2.4);
 107     else
 108         fg = fg / 12.92;
 109
 110     if (fb > 0.04045)
 111         fb = (float) pow((fb + 0.055) / 1.055, 2.4);
 112     else
 113         fb = fb / 12.92;
 114
 115     fr = fr * 100.0;
 116     fg = fg * 100.0;
 117     fb = fb * 100.0;
 118
 119     // Use white = D65
 120     float x = fr * 0.4124 + fg * 0.3576 + fb * 0.1805;
 121     float y = fr * 0.2126 + fg * 0.7152 + fb * 0.0722;
 122     float z = fr * 0.0193 + fg * 0.1192 + fb * 0.9505;
 123
 124     float vx = x /  95.047;
 125     float vy = y / 100.000;
 126     float vz = z / 108.883;
 127
 128     if (vx > 0.008856)
 129         vx = (float) pow(vx, 0.3333);
 130     else
 131         vx = (7.787 * vx) + (16.0 / 116.0);
 132
 133     if (vy > 0.008856)
 134         vy = (float) pow(vy, 0.3333);
 135     else
 136         vy = (7.787 * vy) + (16.0 / 116.0);
 137
 138     if (vz > 0.008856)
 139         vz = (float) pow(vz, 0.3333);
 140     else
 141         vz = (7.787 * vz) + (16.0 / 116.0);
 142
 143     C = 0;
 144     L = 116.0 * vy - 16.0;
 145     A = 500.0 * (vx - vy);
 146     B = 200.0 * (vy - vz);
 147
 148     // Cache for next time
 149     clabLookupTable[rgb] = *this;
 150
 151 }
 152
 153
 154
 155 /**
 156  * Return this CLAB's value a a packed-pixel ARGB value
 157  */
 158 unsigned long CLAB::toRGB()
 159 {
 160     float vy = (L + 16.0) / 116.0;
 161     float vx = A / 500.0 + vy;
 162     float vz = vy - B / 200.0;
 163
 164     float vx3 = vx * vx * vx;
 165     float vy3 = vy * vy * vy;
 166     float vz3 = vz * vz * vz;
 167
 168     if (vy3 > 0.008856)
 169         vy = vy3;
 170     else
 171         vy = (vy - 16.0 / 116.0) / 7.787;
 172
 173     if (vx3 > 0.008856)
 174         vx = vx3;
 175     else
 176         vx = (vx - 16.0 / 116.0) / 7.787;
 177
 178     if (vz3 > 0.008856)
 179         vz = vz3;
 180     else
 181         vz = (vz - 16.0 / 116.0) / 7.787;
 182
 183     float x =  95.047 * vx; //use white = D65
 184     float y = 100.000 * vy;
 185     float z = 108.883 * vz;
 186
 187     vx = x / 100.0;
 188     vy = y / 100.0;
 189     vz = z / 100.0;
 190
 191     float vr =(float)(vx *  3.2406 + vy * -1.5372 + vz * -0.4986);
 192     float vg =(float)(vx * -0.9689 + vy *  1.8758 + vz *  0.0415);
 193     float vb =(float)(vx *  0.0557 + vy * -0.2040 + vz *  1.0570);
 194
 195     if (vr > 0.0031308)
 196         vr = (float)(1.055 * pow(vr, (1.0 / 2.4)) - 0.055);
 197     else
 198         vr = 12.92 * vr;
 199
 200     if (vg > 0.0031308)
 201         vg = (float)(1.055 * pow(vg, (1.0 / 2.4)) - 0.055);
 202     else
 203         vg = 12.92 * vg;
 204
 205     if (vb > 0.0031308)
 206         vb = (float)(1.055 * pow(vb, (1.0 / 2.4)) - 0.055);
 207     else
 208         vb = 12.92 * vb;
 209
 210     return getRGB(0.0, vr, vg, vb);
 211 }
 212
 213
 214 /**
 215  * Computes squared euclidian distance in CLAB space for two colors
 216  * given as RGB values.
 217  */
 218 float CLAB::diffSq(unsigned int rgb1, unsigned int rgb2)
 219 {
 220     CLAB c1(rgb1);
 221     CLAB c2(rgb2);
 222     float euclid=0.0f;
 223     euclid += (c1.L - c2.L) * (c1.L - c2.L);
 224     euclid += (c1.A - c2.A) * (c1.A - c2.A);
 225     euclid += (c1.B - c2.B) * (c1.B - c2.B);
 226     return euclid;
 227 }
 228
 229
 230 /**
 231  * Computes squared euclidian distance in CLAB space for two colors
 232  * given as RGB values.
 233  */
 234 float CLAB::diff(unsigned int rgb0, unsigned int rgb1)
 235 {
 236     return (float) sqrt(diffSq(rgb0, rgb1));
 237 }
 238
 239
 240 //########################################################################
 241 //#  T U P E L
 242 //########################################################################
 243
 244 /**
 245  * Helper class for storing the minimum distances to a cluster centroid
 246  * in background and foreground and the index to the centroids in each
 247  * signature for a given color.
 248  */
 249 class Tupel {
 250 public:
 251
 252     Tupel()
 253         {
 254         minBgDist  = 0.0f;
 255         indexMinBg = 0;
 256         minFgDist  = 0.0f;
 257         indexMinFg = 0;
 258         }
 259     Tupel(float minBgDistArg, long indexMinBgArg,
 260           float minFgDistArg, long indexMinFgArg)
 261         {
 262         minBgDist  = minBgDistArg;
 263         indexMinBg = indexMinBgArg;
 264         minFgDist  = minFgDistArg;
 265         indexMinFg = indexMinFgArg;
 266         }
 267     Tupel(const Tupel &other)
 268         {
 269         minBgDist  = other.minBgDist;
 270         indexMinBg = other.indexMinBg;
 271         minFgDist  = other.minFgDist;
 272         indexMinFg = other.indexMinFg;
 273         }
 274     Tupel &operator=(const Tupel &other)
 275         {
 276         minBgDist  = other.minBgDist;
 277         indexMinBg = other.indexMinBg;
 278         minFgDist  = other.minFgDist;
 279         indexMinFg = other.indexMinFg;
 280         return *this;
 281         }
 282     virtual ~Tupel()
 283         {}
 284
 285     float minBgDist;
 286     long  indexMinBg;
 287     float minFgDist;
 288     long  indexMinFg;
 289
 290  };
 291
 292
 293
 294 //########################################################################
 295 //#  S I O X    I M A G E
 296 //########################################################################
 297
 298 /**
 299  *  Create an image with the given width and height
 300  */
 301 SioxImage::SioxImage(unsigned int widthArg, unsigned int heightArg)
 302 {
 303     init(width, height);
 304 }
 305
 306 /**
 307  *  Copy constructor
 308  */
 309 SioxImage::SioxImage(const SioxImage &other)
 310 {
 311     pixdata = NULL;
 312     cmdata  = NULL;
 313     assign(other);
 314 }
 315
 316 /**
 317  *  Assignment
 318  */
 319 SioxImage &SioxImage::operator=(const SioxImage &other)
 320 {
 321     assign(other);
 322     return *this;
 323 }
 324
 325
 326 /**
 327  * Clean up after use.
 328  */
 329 SioxImage::~SioxImage()
 330 {
 331     if (pixdata) delete[] pixdata;
 332     if (cmdata)  delete[] cmdata;
 333 }
 334
 335 /**
 336  * Returns true if the previous operation on this image
 337  * was successful, else false.
 338  */
 339 bool SioxImage::isValid()
 340 {
 341     return valid;
 342 }
 343
 344 /**
 345  * Sets whether an operation was successful, and whether
 346  * this image should be considered a valid one.
 347  * was successful, else false.
 348  */
 349 void SioxImage::setValid(bool val)
 350 {
 351     valid = val;
 352 }
 353
 354
 355 /**
 356  * Set a pixel at the x,y coordinates to the given value.
 357  * If the coordinates are out of range, do nothing.
 358  */
 359 void SioxImage::setPixel(unsigned int x,
 360                          unsigned int y,
 361                          unsigned int pixval)
 362 {
 363     if (x > width || y > height)
 364         return;
 365     unsigned long offset = width * y + x;
 366     pixdata[offset] = pixval;
 367 }
 368
 369 /**
 370  * Set a pixel at the x,y coordinates to the given r, g, b values.
 371  * If the coordinates are out of range, do nothing.
 372  */
 373 void SioxImage::setPixel(unsigned int x, unsigned int y,
 374                          unsigned int a,
 375                          unsigned int r,
 376                          unsigned int g,
 377                          unsigned int b)
 378 {
 379     if (x > width || y > height)
 380         return;
 381     unsigned long offset = width * y + x;
 382     unsigned int pixval = ((a << 24) & 0xff000000) |
 383                           ((r << 16) & 0x00ff0000) |
 384                           ((g <<  8) & 0x0000ff00) |
 385                           ((b      ) & 0x000000ff);
 386     pixdata[offset] = pixval;
 387 }
 388
 389
 390
 391 /**
 392  *  Get a pixel at the x,y coordinates given.  If
 393  *  the coordinates are out of range, return 0;
 394  */
 395 unsigned int SioxImage::getPixel(unsigned int x, unsigned int y)
 396 {
 397     if (x > width || y > height)
 398         return 0L;
 399     unsigned long offset = width * y + x;
 400     return pixdata[offset];
 401 }
 402
 403 /**
 404  *  Return the image data buffer
 405  */
 406 unsigned int *SioxImage::getImageData()
 407 {
 408     return pixdata;
 409 }
 410
 411 /**
 412  * Set a confidence value at the x,y coordinates to the given value.
 413  * If the coordinates are out of range, do nothing.
 414  */
 415 void SioxImage::setConfidence(unsigned int x,
 416                               unsigned int y,
 417                               float confval)
 418 {
 419     if (x > width || y > height)
 420         return;
 421     unsigned long offset = width * y + x;
 422     cmdata[offset] = confval;
 423 }
 424
 425 /**
 426  *  Get a confidence valueat the x,y coordinates given.  If
 427  *  the coordinates are out of range, return 0;
 428  */
 429 float SioxImage::getConfidence(unsigned int x, unsigned int y)
 430 {
 431     if (x > width || y > height)
 432         return 0.0;
 433     unsigned long offset = width * y + x;
 434     return cmdata[offset];
 435 }
 436
 437 /**
 438  *  Return the confidence data buffer
 439  */
 440 float *SioxImage::getConfidenceData()
 441 {
 442     return cmdata;
 443 }
 444
 445
 446 /**
 447  * Return the width of this image
 448  */
 449 int SioxImage::getWidth()
 450 {
 451     return width;
 452 }
 453
 454 /**
 455  * Return the height of this image
 456  */
 457 int SioxImage::getHeight()
 458 {
 459     return height;
 460 }
 461
 462 /**
 463  * Initialize values.  Used by constructors
 464  */
 465 void SioxImage::init(unsigned int widthArg, unsigned int heightArg)
 466 {
 467     valid     = true;
 468     width     = widthArg;
 469     height    = heightArg;
 470     imageSize = width * height;
 471     pixdata   = new unsigned int[imageSize];
 472     cmdata    = new float[imageSize];
 473     for (unsigned long i=0 ; i<imageSize ; i++)
 474         {
 475         pixdata[i] = 0;
 476         cmdata[i]  = 0.0;
 477         }
 478 }
 479
 480 /**
 481  * Assign values to that of another
 482  */
 483 void SioxImage::assign(const SioxImage &other)
 484 {
 485     if (pixdata) delete[] pixdata;
 486     if (cmdata)  delete[] cmdata;
 487     valid     = other.valid;
 488     width     = other.width;
 489     height    = other.height;
 490     imageSize = width * height;
 491     pixdata   = new unsigned int[imageSize];
 492     cmdata    = new float[imageSize];
 493     for (unsigned long i=0 ; i<imageSize ; i++)
 494         {
 495         pixdata[i] = other.pixdata[i];
 496         cmdata[i]  = other.cmdata[i];
 497         }
 498 }
 499
 500
 501 /**
 502  * Write the image to a PPM file
 503  */
 504 bool SioxImage::writePPM(const std::string fileName)
 505 {
 506
 507     FILE *f = fopen(fileName.c_str(), "wb");
 508     if (!f)
 509         return false;
 510
 511     fprintf(f, "P6 %d %d 255\n", width, height);
 512
 513     for (unsigned int y=0 ; y<height; y++)
 514         {
 515         for (unsigned int x=0 ; x<width ; x++)
 516             {
 517             unsigned int rgb = getPixel(x, y);
 518             //unsigned int alpha = (rgb>>24) & 0xff;
 519             unsigned int r = ((rgb>>16) & 0xff);
 520             unsigned int g = ((rgb>> 8) & 0xff);
 521             unsigned int b = ((rgb    ) & 0xff);
 522             fputc((unsigned char) r, f);
 523             fputc((unsigned char) g, f);
 524             fputc((unsigned char) b, f);
 525             }
 526         }
 527     fclose(f);
 528     return true;
 529 }
 530
 531
 532 #ifdef HAVE_GLIB
 533
 534 /**
 535  * Create an image from a GdkPixbuf
 536  */
 537 SioxImage::SioxImage(GdkPixbuf *buf)
 538 {
 539     if (!buf)
 540         return;
 541
 542     unsigned int width  = gdk_pixbuf_get_width(buf);
 543     unsigned int height = gdk_pixbuf_get_height(buf);
 544     init(width, height); //DO THIS NOW!!
 545
 546
 547     guchar *pixldata    = gdk_pixbuf_get_pixels(buf);
 548     int rowstride       = gdk_pixbuf_get_rowstride(buf);
 549     int n_channels      = gdk_pixbuf_get_n_channels(buf);
 550
 551     //### Fill in the cells with RGB values
 552     int row  = 0;
 553     for (unsigned int y=0 ; y<height ; y++)
 554         {
 555         guchar *p = pixldata + row;
 556         for (unsigned int x=0 ; x<width ; x++)
 557             {
 558             int r     = (int)p[0];
 559             int g     = (int)p[1];
 560             int b     = (int)p[2];
 561             int alpha = (int)p[3];
 562
 563             setPixel(x, y, alpha, r, g, b);
 564             p += n_channels;
 565             }
 566         row += rowstride;
 567         }
 568
 569 }
 570
 571
 572 /**
 573  * Create a GdkPixbuf from this image
 574  */
 575 GdkPixbuf *SioxImage::getGdkPixbuf()
 576 {
 577     guchar *pixdata = (guchar *)
 578           malloc(sizeof(guchar) * width * height * 3);
 579     if (!pixdata)
 580         return NULL;
 581
 582     int n_channels = 3;
 583     int rowstride  = width * 3;
 584
 585     GdkPixbuf *buf = gdk_pixbuf_new_from_data(pixdata, GDK_COLORSPACE_RGB,
 586                         0, 8, width, height,
 587                         rowstride, NULL, NULL);
 588
 589     //### Fill in the cells with RGB values
 590     int row  = 0;
 591     for (unsigned int y=0 ; y < height ; y++)
 592         {
 593         guchar *p = pixdata + row;
 594         for (unsigned x=0 ; x < width ; x++)
 595             {
 596             unsigned int rgb = getPixel(x, y);
 597             p[0] = (rgb >> 16) & 0xff;
 598             p[1] = (rgb >>  8) & 0xff;
 599             p[2] = (rgb      ) & 0xff;
 600             p += n_channels;
 601             }
 602         row += rowstride;
 603         }
 604
 605     return buf;
 606 }
 607
 608 #endif /* GLIB */
 609
 610
 611
 612
 613 //########################################################################
 614 //#  S I O X
 615 //########################################################################
 616
 617 //##############
 618 //## PUBLIC
 619 //##############
 620
 621 /**
 622  * Confidence corresponding to a certain foreground region (equals one).
 623  */
 624 const float Siox::CERTAIN_FOREGROUND_CONFIDENCE=1.0f;
 625
 626 /**
 627  * Confidence for a region likely being foreground.
 628  */
 629 const float Siox::FOREGROUND_CONFIDENCE=0.8f;
 630
 631 /**
 632  * Confidence for foreground or background type being equally likely.
 633  */
 634 const float Siox::UNKNOWN_REGION_CONFIDENCE=0.5f;
 635
 636 /**
 637  * Confidence for a region likely being background.
 638  */
 639 const float Siox::BACKGROUND_CONFIDENCE=0.1f;
 640
 641 /**
 642  * Confidence corresponding to a certain background reagion (equals zero).
 643  */
 644 const float Siox::CERTAIN_BACKGROUND_CONFIDENCE=0.0f;
 645
 646 /**
 647  *  Construct a Siox engine
 648  */
 649 Siox::Siox()
 650 {
 651     init();
 652 }
 653
 654
 655 /**
 656  *
 657  */
 658 Siox::~Siox()
 659 {
 660     cleanup();
 661 }
 662
 663
 664 /**
 665  * Error logging
 666  */
 667 void Siox::error(char *fmt, ...)
 668 {
 669     char msgbuf[256];
 670     va_list args;
 671     va_start(args, fmt);
 672     vsnprintf(msgbuf, 255, fmt, args);
 673     va_end(args) ;
 674 #ifdef HAVE_GLIB
 675     g_warning("Siox error: %s\n", msgbuf);
 676 #else
 677     fprintf(stderr, "Siox error: %s\n", msgbuf);
 678 #endif
 679 }
 680
 681 /**
 682  * Trace logging
 683  */
 684 void Siox::trace(char *fmt, ...)
 685 {
 686     char msgbuf[256];
 687     va_list args;
 688     va_start(args, fmt);
 689     vsnprintf(msgbuf, 255, fmt, args);
 690     va_end(args) ;
 691 #ifdef HAVE_GLIB
 692     g_message("Siox: %s\n", msgbuf);
 693 #else
 694     fprintf(stdout, "Siox: %s\n", msgbuf);
 695 #endif
 696 }
 697
 698
 699
 700
 701 /**
 702  *  Extract the foreground of the original image, according
 703  *  to the values in the confidence matrix.
 704  */
 705 SioxImage Siox::extractForeground(const SioxImage &originalImage,
 706                                   unsigned int backgroundFillColor)
 707 {
 708     init();
 709
 710     SioxImage workImage = originalImage;
 711
 712     //# fetch some info from the image
 713     width      = workImage.getWidth();
 714     height     = workImage.getHeight();
 715     pixelCount = width * height;
 716     image      = workImage.getImageData();
 717     cm         = workImage.getConfidenceData();
 718     labelField = new int[pixelCount];
 719
 720     trace("### Creating signatures");
 721
 722     //#### create color signatures
 723     std::vector<CLAB> knownBg;
 724     std::vector<CLAB> knownFg;
 725     for (int x = 0 ; x < workImage.getWidth() ; x++)
 726         for (int y = 0 ; y < workImage.getHeight() ; y++)
 727             {
 728             float cm = workImage.getConfidence(x, y);
 729             unsigned int pix = workImage.getPixel(x, y);
 730             if (cm <= BACKGROUND_CONFIDENCE)
 731                 knownBg.push_back(pix); //note: uses CLAB(rgb)
 732             else if (cm >= FOREGROUND_CONFIDENCE)
 733                 knownFg.push_back(pix);
 734             }
 735
 736     trace("knownBg:%d knownFg:%d", knownBg.size(), knownFg.size());
 737
 738
 739     std::vector<CLAB> bgSignature ;
 740     if (!colorSignature(knownBg, bgSignature, 3))
 741         {
 742         error("Could not create background signature");
 743         workImage.setValid(false);
 744         return workImage;
 745         }
 746     std::vector<CLAB> fgSignature ;
 747     if (!colorSignature(knownFg, fgSignature, 3))
 748         {
 749         error("Could not create foreground signature");
 750         delete[] labelField;
 751         workImage.setValid(false);
 752         return workImage;
 753         }
 754
 755     //trace("### bgSignature:%d", bgSignature.size());
 756
 757     if (bgSignature.size() < 1)
 758         {
 759         // segmentation impossible
 760         error("Signature size is < 1.  Segmentation is impossible");
 761         delete[] labelField;
 762         workImage.setValid(false);
 763         return workImage;
 764         }
 765
 766     // classify using color signatures,
 767     // classification cached in hashmap for drb and speedup purposes
 768     std::map<unsigned int, Tupel> hs;
 769
 770     for (unsigned int i=0; i<pixelCount; i++)
 771         {
 772         if (cm[i] >= FOREGROUND_CONFIDENCE)
 773             {
 774             cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
 775             }
 776         else if (cm[i] <= BACKGROUND_CONFIDENCE)
 777             {
 778             cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
 779             }
 780         else // somewhere in between
 781             {
 782             bool isBackground = true;
 783             std::map<unsigned int, Tupel>::iterator iter = hs.find(i);
 784             if (iter != hs.end()) //found
 785                 {
 786                 Tupel tupel = iter->second;
 787                 isBackground = tupel.minBgDist <= tupel.minFgDist;
 788                 }
 789             else
 790                 {
 791                 CLAB lab(image[i]);
 792                 float minBg = sqrEuclidianDist(lab, bgSignature[0]);
 793                 int minIndex=0;
 794                 for (unsigned int j=1; j<bgSignature.size() ; j++)
 795                     {
 796                     float d = sqrEuclidianDist(lab, bgSignature[j]);
 797                     if (d<minBg)
 798                         {
 799                         minBg    = d;
 800                         minIndex = j;
 801                         }
 802                     }
 803                 Tupel tupel(0.0f, 0,  0.0f, 0);
 804                 tupel.minBgDist  = minBg;
 805                 tupel.indexMinBg = minIndex;
 806                 float minFg = 1.0e6f;
 807                 minIndex = -1;
 808                 for (unsigned int j = 0 ; j < fgSignature.size() ; j++)
 809                     {
 810                     float d = sqrEuclidianDist(lab, fgSignature[j]);
 811                     if (d < minFg)
 812                         {
 813                         minFg    = d;
 814                         minIndex = j;
 815                         }
 816                     }
 817                 tupel.minFgDist  = minFg;
 818                 tupel.indexMinFg = minIndex;
 819                 if (fgSignature.size() == 0)
 820                     {
 821                     isBackground=(minBg <= clusterSize);
 822                     // remove next line to force behaviour of old algorithm
 823                     //error("foreground signature does not exist");
 824                     //delete[] labelField;
 825                     //workImage.setValid(false);
 826                     //return workImage;
 827                     }
 828                 else
 829                     {
 830                     isBackground = minBg < minFg;
 831                     }
 832                 hs[image[i]] = tupel;
 833                 }
 834
 835             if (isBackground)
 836                 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
 837             else
 838                 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
 839             }
 840         }
 841
 842     trace("### postProcessing");
 843
 844     //## postprocessing
 845     smooth(cm, width, height, 0.33f, 0.33f, 0.33f); // average
 846     normalizeMatrix(cm, pixelCount);
 847     erode(cm, width, height);
 848     keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.0/*sizeFactorToKeep*/);
 849
 850     for (int i=0; i < 2/*smoothness*/; i++)
 851         smooth(cm, width, height, 0.33f, 0.33f, 0.33f); // average
 852
 853     normalizeMatrix(cm, pixelCount);
 854
 855     for (unsigned int i=0; i < pixelCount; i++)
 856         {
 857         if (cm[i] >= UNKNOWN_REGION_CONFIDENCE)
 858             cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
 859         else
 860             cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
 861         }
 862
 863     keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.5/*sizeFactorToKeep*/);
 864     fillColorRegions();
 865     dilate(cm, width, height);
 866
 867     delete[] labelField;
 868
 869     //#### Yaay.  We are done.  Now clear everything but the background
 870     for (unsigned int y = 0 ; y < height ; y++)
 871         for (unsigned int x = 0 ; x < width ; x++)
 872             {
 873             float conf = workImage.getConfidence(x, y);
 874             if (conf < FOREGROUND_CONFIDENCE)
 875                 {
 876                 workImage.setPixel(x, y, backgroundFillColor);
 877                 }
 878             }
 879
 880     trace("### Done");
 881     return workImage;
 882 }
 883
 884
 885
 886 //##############
 887 //## PRIVATE
 888 //##############
 889
 890 /**
 891  *  Initialize the Siox engine to its 'pristine' state.
 892  *  Performed at the beginning of extractForeground().
 893  */
 894 void Siox::init()
 895 {
 896     limits[0] = 0.64f;
 897     limits[1] = 1.28f;
 898     limits[2] = 2.56f;
 899
 900     float negLimits[3];
 901     negLimits[0] = -limits[0];
 902     negLimits[1] = -limits[1];
 903     negLimits[2] = -limits[2];
 904
 905     clusterSize = sqrEuclidianDist(limits, 3, negLimits);
 906 }
 907
 908
 909 /**
 910  *  Clean up any debris from processing.
 911  */
 912 void Siox::cleanup()
 913 {
 914 }
 915
 916
 917
 918
 919 /**
 920  *  Stage 1 of the color signature work.  'dims' will be either
 921  *  2 for grays, or 3 for colors
 922  */
 923 void Siox::colorSignatureStage1(CLAB *points,
 924                                 unsigned int leftBase,
 925                                 unsigned int rightBase,
 926                                 unsigned int recursionDepth,
 927                                 unsigned int *clusterCount,
 928                                 const unsigned int dims)
 929 {
 930
 931     unsigned int currentDim = recursionDepth % dims;
 932     CLAB point = points[leftBase];
 933     float min = point(currentDim);
 934     float max = min;
 935
 936     for (unsigned int i = leftBase + 1; i < rightBase ; i++)
 937         {
 938         point = points[i];
 939         float curval = point(currentDim);
 940         if (curval < min) min = curval;
 941         if (curval > max) max = curval;
 942         }
 943
 944     //Do the Rubner-rule split (sounds like a dance)
 945     if (max - min > limits[currentDim])
 946         {
 947         float pivotPoint = (min + max) / 2.0; //average
 948         unsigned int left  = leftBase;
 949         unsigned int right = rightBase - 1;
 950
 951         //# partition points according to the dimension
 952         while (true)
 953             {
 954             while ( true )
 955                 {
 956                 point = points[left];
 957                 if (point(currentDim) > pivotPoint)
 958                     break;
 959                 left++;
 960                 }
 961             while ( true )
 962                 {
 963                 point = points[right];
 964                 if (point(currentDim) <= pivotPoint)
 965                     break;
 966                 right--;
 967                 }
 968
 969             if (left > right)
 970                 break;
 971
 972             point = points[left];
 973             points[left] = points[right];
 974             points[right] = point;
 975
 976             left++;
 977             right--;
 978             }
 979
 980         //# Recurse and create sub-trees
 981         colorSignatureStage1(points, leftBase, left,
 982                  recursionDepth + 1, clusterCount, dims);
 983         colorSignatureStage1(points, left, rightBase,
 984                  recursionDepth + 1, clusterCount, dims);
 985         }
 986     else
 987         {
 988         //create a leaf
 989         CLAB newpoint;
 990
 991         newpoint.C = rightBase - leftBase;
 992
 993         for (; leftBase < rightBase ; leftBase++)
 994             {
 995             newpoint.add(points[leftBase]);
 996             }
 997
 998         //printf("clusters:%d\n", *clusters);
 999
1000         if (newpoint.C != 0)
1001             newpoint.mul(1.0 / (float)newpoint.C);
1002         points[*clusterCount] = newpoint;
1003         (*clusterCount)++;
1004         }
1005 }
1006
1007
1008
1009 /**
1010  *  Stage 2 of the color signature work
1011  */
1012 void Siox::colorSignatureStage2(CLAB         *points,
1013                                 unsigned int leftBase,
1014                                 unsigned int rightBase,
1015                                 unsigned int recursionDepth,
1016                                 unsigned int *clusterCount,
1017                                 const float  threshold,
1018                                 const unsigned int dims)
1019 {
1020
1021
1022     unsigned int currentDim = recursionDepth % dims;
1023     CLAB point = points[leftBase];
1024     float min = point(currentDim);
1025     float max = min;
1026
1027     for (unsigned int i = leftBase+ 1; i < rightBase; i++)
1028         {
1029         point = points[i];
1030         float curval = point(currentDim);
1031         if (curval < min) min = curval;
1032         if (curval > max) max = curval;
1033         }
1034
1035     //Do the Rubner-rule split (sounds like a dance)
1036     if (max - min > limits[currentDim])
1037         {
1038         float pivotPoint = (min + max) / 2.0; //average
1039         unsigned int left  = leftBase;
1040         unsigned int right = rightBase - 1;
1041
1042         //# partition points according to the dimension
1043         while (true)
1044             {
1045             while ( true )
1046                 {
1047                 point = points[left];
1048                 if (point(currentDim) > pivotPoint)
1049                     break;
1050                 left++;
1051                 }
1052             while ( true )
1053                 {
1054                 point = points[right];
1055                 if (point(currentDim) <= pivotPoint)
1056                     break;
1057                 right--;
1058                 }
1059
1060             if (left > right)
1061                 break;
1062
1063             point = points[left];
1064             points[left] = points[right];
1065             points[right] = point;
1066
1067             left++;
1068             right--;
1069             }
1070
1071         //# Recurse and create sub-trees
1072         colorSignatureStage2(points, leftBase, left,
1073                  recursionDepth + 1, clusterCount, threshold, dims);
1074         colorSignatureStage2(points, left, rightBase,
1075                  recursionDepth + 1, clusterCount, threshold, dims);
1076         }
1077     else
1078         {
1079         //### Create a leaf
1080         unsigned int sum = 0;
1081         for (unsigned int i = leftBase; i < rightBase; i++)
1082             sum += points[i].C;
1083
1084         if ((float)sum >= threshold)
1085             {
1086             float scale = (float)(rightBase - leftBase);
1087             CLAB newpoint;
1088
1089             for (; leftBase < rightBase; leftBase++)
1090                 newpoint.add(points[leftBase]);
1091
1092             if (scale != 0.0)
1093                 newpoint.mul(1.0 / scale);
1094             points[*clusterCount] = newpoint;
1095             (*clusterCount)++;
1096             }
1097       }
1098 }
1099
1100
1101
1102 /**
1103  *  Main color signature method
1104  */
1105 bool Siox::colorSignature(const std::vector<CLAB> &inputVec,
1106                           std::vector<CLAB> &result,
1107                           const unsigned int dims)
1108 {
1109
1110     unsigned int length = inputVec.size();
1111
1112     if (length < 1) // no error. just don't do anything
1113         return true;
1114
1115     CLAB *input = (CLAB *) malloc(length * sizeof(CLAB));
1116
1117     if (!input)
1118         {
1119         error("Could not allocate buffer for signature");
1120         return false;
1121         }
1122     for (unsigned int i=0 ; i < length ; i++)
1123         input[i] = inputVec[i];
1124
1125     unsigned int stage1length = 0;
1126     colorSignatureStage1(input, 0, length, 0,
1127                    &stage1length, dims);
1128
1129     unsigned int stage2length = 0;
1130     colorSignatureStage2(input, 0, stage1length, 0,
1131                   &stage2length, length * 0.001, dims);
1132
1133     result.clear();
1134     for (unsigned int i=0 ; i < stage2length ; i++)
1135         result.push_back(input[i]);
1136
1137     free(input);
1138
1139     return true;
1140 }
1141
1142
1143
1144 /**
1145  *
1146  */
1147 void Siox::keepOnlyLargeComponents(float threshold,
1148                                    double sizeFactorToKeep)
1149 {
1150     for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1151         labelField[idx] = -1;
1152
1153     int curlabel = 0;
1154     int maxregion= 0;
1155     int maxblob  = 0;
1156
1157     // slow but easy to understand:
1158     std::vector<int> labelSizes;
1159     for (unsigned long i=0 ; i<pixelCount ; i++)
1160         {
1161         int regionCount = 0;
1162         if (labelField[i] == -1 && cm[i] >= threshold)
1163             {
1164             regionCount = depthFirstSearch(i, threshold, curlabel++);
1165             labelSizes.push_back(regionCount);
1166             }
1167
1168         if (regionCount>maxregion)
1169             {
1170             maxregion = regionCount;
1171             maxblob   = curlabel-1;
1172             }
1173         }
1174
1175     for (unsigned int i=0 ; i<pixelCount ; i++)
1176         {
1177         if (labelField[i] != -1)
1178             {
1179             // remove if the component is to small
1180             if (labelSizes[labelField[i]] * sizeFactorToKeep < maxregion)
1181                 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
1182
1183             // add maxblob always to foreground
1184             if (labelField[i] == maxblob)
1185                 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1186             }
1187         }
1188
1189 }
1190
1191
1192
1193 int Siox::depthFirstSearch(int startPos,
1194               float threshold, int curLabel)
1195 {
1196     // stores positions of labeled pixels, where the neighbours
1197     // should still be checked for processing:
1198
1199     //trace("startPos:%d threshold:%f curLabel:%d",
1200     //     startPos, threshold, curLabel);
1201
1202     std::vector<int> pixelsToVisit;
1203     int componentSize = 0;
1204
1205     if (labelField[startPos]==-1 && cm[startPos]>=threshold)
1206         {
1207         labelField[startPos] = curLabel;
1208         componentSize++;
1209         pixelsToVisit.push_back(startPos);
1210         }
1211
1212
1213     while (pixelsToVisit.size() > 0)
1214         {
1215         int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1216         pixelsToVisit.erase(pixelsToVisit.end() - 1);
1217         unsigned int x = pos % width;
1218         unsigned int y = pos / width;
1219
1220         // check all four neighbours
1221         int left = pos - 1;
1222         if (((int)x)-1>=0 && labelField[left]==-1 && cm[left]>=threshold)
1223             {
1224             labelField[left]=curLabel;
1225             componentSize++;
1226             pixelsToVisit.push_back(left);
1227             }
1228
1229         int right = pos + 1;
1230         if (x+1 < width && labelField[right]==-1 && cm[right]>=threshold)
1231             {
1232             labelField[right]=curLabel;
1233             componentSize++;
1234             pixelsToVisit.push_back(right);
1235             }
1236
1237         int top = pos - width;
1238         if (((int)y)-1>=0 && labelField[top]==-1 && cm[top]>=threshold)
1239             {
1240             labelField[top]=curLabel;
1241             componentSize++;
1242             pixelsToVisit.push_back(top);
1243             }
1244
1245         int bottom = pos + width;
1246         if (y+1 < height && labelField[bottom]==-1
1247                          && cm[bottom]>=threshold)
1248             {
1249             labelField[bottom]=curLabel;
1250             componentSize++;
1251             pixelsToVisit.push_back(bottom);
1252             }
1253
1254         }
1255     return componentSize;
1256 }
1257
1258
1259
1260 /**
1261  *
1262  */
1263 void Siox::fillColorRegions()
1264 {
1265     for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1266         labelField[idx] = -1;
1267
1268     //int maxRegion=0; // unused now
1269     std::vector<int> pixelsToVisit;
1270     for (unsigned long i=0; i<pixelCount; i++)
1271         { // for all pixels
1272         if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE)
1273             {
1274             continue; // already visited or bg
1275             }
1276
1277         unsigned int  origColor = image[i];
1278         unsigned long curLabel  = i+1;
1279         labelField[i]           = curLabel;
1280         cm[i]                   = CERTAIN_FOREGROUND_CONFIDENCE;
1281
1282         // int componentSize = 1;
1283         pixelsToVisit.push_back(i);
1284         // depth first search to fill region
1285         while (pixelsToVisit.size() > 0)
1286             {
1287             int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1288             pixelsToVisit.erase(pixelsToVisit.end() - 1);
1289             unsigned int x=pos % width;
1290             unsigned int y=pos / width;
1291             // check all four neighbours
1292             int left = pos-1;
1293             if (((int)x)-1 >= 0 && labelField[left] == -1
1294                         && CLAB::diff(image[left], origColor)<1.0)
1295                 {
1296                 labelField[left]=curLabel;
1297                 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1298                 // ++componentSize;
1299                 pixelsToVisit.push_back(left);
1300                 }
1301             int right = pos+1;
1302             if (x+1 < width && labelField[right]==-1
1303                         && CLAB::diff(image[right], origColor)<1.0)
1304                 {
1305                 labelField[right]=curLabel;
1306                 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1307                 // ++componentSize;
1308                 pixelsToVisit.push_back(right);
1309                 }
1310             int top = pos - width;
1311             if (((int)y)-1>=0 && labelField[top]==-1
1312                         && CLAB::diff(image[top], origColor)<1.0)
1313                 {
1314                 labelField[top]=curLabel;
1315                 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1316                 // ++componentSize;
1317                 pixelsToVisit.push_back(top);
1318                 }
1319             int bottom = pos + width;
1320             if (y+1 < height && labelField[bottom]==-1
1321                         && CLAB::diff(image[bottom], origColor)<1.0)
1322                 {
1323                 labelField[bottom]=curLabel;
1324                 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1325                 // ++componentSize;
1326                 pixelsToVisit.push_back(bottom);
1327                 }
1328             }
1329         //if (componentSize>maxRegion) {
1330         //    maxRegion=componentSize;
1331         //}
1332         }
1333 }
1334
1335
1336
1337
1338 /**
1339  * Applies the morphological dilate operator.
1340  *
1341  * Can be used to close small holes in the given confidence matrix.
1342  */
1343 void Siox::dilate(float *cm, int xres, int yres)
1344 {
1345
1346     for (int y=0; y<yres; y++)
1347         {
1348         for (int x=0; x<xres-1; x++)
1349              {
1350              int idx=(y*xres)+x;
1351              if (cm[idx+1]>cm[idx])
1352                  cm[idx]=cm[idx+1];
1353              }
1354         }
1355
1356     for (int y=0; y<yres; y++)
1357         {
1358         for (int x=xres-1; x>=1; x--)
1359             {
1360             int idx=(y*xres)+x;
1361             if (cm[idx-1]>cm[idx])
1362                 cm[idx]=cm[idx-1];
1363             }
1364         }
1365
1366     for (int y=0; y<yres-1; y++)
1367         {
1368         for (int x=0; x<xres; x++)
1369             {
1370             int idx=(y*xres)+x;
1371             if (cm[((y+1)*xres)+x] > cm[idx])
1372                 cm[idx]=cm[((y+1)*xres)+x];
1373             }
1374         }
1375
1376     for (int y=yres-1; y>=1; y--)
1377         {
1378         for (int x=0; x<xres; x++)
1379             {
1380             int idx=(y*xres)+x;
1381             if (cm[((y-1)*xres)+x] > cm[idx])
1382                 cm[idx]=cm[((y-1)*xres)+x];
1383             }
1384         }
1385 }
1386
1387 /**
1388  * Applies the morphological erode operator.
1389  */
1390 void Siox::erode(float *cm, int xres, int yres)
1391 {
1392     for (int y=0; y<yres; y++)
1393         {
1394         for (int x=0; x<xres-1; x++)
1395             {
1396             int idx=(y*xres)+x;
1397             if (cm[idx+1] < cm[idx])
1398                 cm[idx]=cm[idx+1];
1399             }
1400         }
1401     for (int y=0; y<yres; y++)
1402         {
1403         for (int x=xres-1; x>=1; x--)
1404             {
1405             int idx=(y*xres)+x;
1406             if (cm[idx-1] < cm[idx])
1407                 cm[idx]=cm[idx-1];
1408             }
1409         }
1410     for (int y=0; y<yres-1; y++)
1411         {
1412         for (int x=0; x<xres; x++)
1413             {
1414             int idx=(y*xres)+x;
1415             if (cm[((y+1)*xres)+x] < cm[idx])
1416                 cm[idx]=cm[((y+1)*xres)+x];
1417             }
1418         }
1419     for (int y=yres-1; y>=1; y--)
1420         {
1421         for (int x=0; x<xres; x++)
1422             {
1423             int idx=(y*xres)+x;
1424             if (cm[((y-1)*xres)+x] < cm[idx])
1425                 cm[idx]=cm[((y-1)*xres)+x];
1426             }
1427         }
1428 }
1429
1430
1431
1432
1433 /**
1434  * Normalizes the matrix to values to [0..1].
1435  */
1436 void Siox::normalizeMatrix(float *cm, int cmSize)
1437 {
1438     float max= -1000000.0f;
1439     for (int i=0; i<cmSize; i++)
1440         if (max<cm[i] > max)
1441             max=cm[i];
1442
1443     if (max<=0.0 || max==1.0)
1444         return;
1445
1446     float alpha=1.00f/max;
1447     premultiplyMatrix(alpha, cm, cmSize);
1448 }
1449
1450 /**
1451  * Multiplies matrix with the given scalar.
1452  */
1453 void Siox::premultiplyMatrix(float alpha, float *cm, int cmSize)
1454 {
1455     for (int i=0; i<cmSize; i++)
1456         cm[i]=alpha*cm[i];
1457 }
1458
1459 /**
1460  * Blurs confidence matrix with a given symmetrically weighted kernel.
1461  *
1462  * In the standard case confidence matrix entries are between 0...1 and
1463  * the weight factors sum up to 1.
1464  */
1465 void Siox::smooth(float *cm, int xres, int yres,
1466                   float f1, float f2, float f3)
1467 {
1468     for (int y=0; y<yres; y++)
1469         {
1470         for (int x=0; x<xres-2; x++)
1471             {
1472             int idx=(y*xres)+x;
1473             cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
1474             }
1475         }
1476     for (int y=0; y<yres; y++)
1477         {
1478         for (int x=xres-1; x>=2; x--)
1479             {
1480             int idx=(y*xres)+x;
1481             cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
1482             }
1483         }
1484     for (int y=0; y<yres-2; y++)
1485         {
1486         for (int x=0; x<xres; x++)
1487             {
1488             int idx=(y*xres)+x;
1489             cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
1490             }
1491         }
1492     for (int y=yres-1; y>=2; y--)
1493         {
1494         for (int x=0; x<xres; x++)
1495             {
1496             int idx=(y*xres)+x;
1497             cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
1498             }
1499         }
1500 }
1501
1502 /**
1503  * Squared Euclidian distance of p and q.
1504  */
1505 float Siox::sqrEuclidianDist(float *p, int pSize, float *q)
1506 {
1507     float sum=0.0;
1508     for (int i=0; i<pSize; i++)
1509         {
1510         float v = p[i] - q[i];
1511         sum += v*v;
1512         }
1513     return sum;
1514 }
1515
1516 /**
1517  * Squared Euclidian distance of p and q.
1518  */
1519 float Siox::sqrEuclidianDist(const CLAB &p, const CLAB &q)
1520 {
1521     float sum=0;
1522     sum += (p.L - q.L) * (p.L - q.L);
1523     sum += (p.A - q.A) * (p.A - q.A);
1524     sum += (p.B - q.B) * (p.B - q.B);
1525     return sum;
1526 }
1527
1528
1529
1530
1531
1532
1533
1534
1535 } // namespace siox
1536 } // namespace org
1537
1538 //########################################################################
1539 //#  E N D    O F    F I L E
1540 //########################################################################
1541