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     bool has_alpha = true;
 578     int n_channels = has_alpha ? 4 : 3;
 579
 580     guchar *pixdata = (guchar *)
 581           malloc(sizeof(guchar) * width * height * n_channels);
 582     if (!pixdata)
 583         return NULL;
 584
 585     int rowstride  = width * n_channels;
 586
 587     GdkPixbuf *buf = gdk_pixbuf_new_from_data(pixdata,
 588                         GDK_COLORSPACE_RGB,
 589                         has_alpha, 8, width, height,
 590                         rowstride, NULL, NULL);
 591
 592     //### Fill in the cells with RGB values
 593     int row  = 0;
 594     for (unsigned int y=0 ; y < height ; y++)
 595         {
 596         guchar *p = pixdata + row;
 597         for (unsigned x=0 ; x < width ; x++)
 598             {
 599             unsigned int rgb = getPixel(x, y);
 600             p[0] = (rgb >> 16) & 0xff;//r
 601             p[1] = (rgb >>  8) & 0xff;//g
 602             p[2] = (rgb      ) & 0xff;//b
 603             if ( n_channels > 3 )
 604                 {
 605                 p[3] = (rgb >> 24) & 0xff;//a
 606                 }
 607             p += n_channels;
 608             }
 609         row += rowstride;
 610         }
 611
 612     return buf;
 613 }
 614
 615 #endif /* GLIB */
 616
 617
 618
 619
 620 //########################################################################
 621 //#  S I O X
 622 //########################################################################
 623
 624 //##############
 625 //## PUBLIC
 626 //##############
 627
 628 /**
 629  * Confidence corresponding to a certain foreground region (equals one).
 630  */
 631 const float Siox::CERTAIN_FOREGROUND_CONFIDENCE=1.0f;
 632
 633 /**
 634  * Confidence for a region likely being foreground.
 635  */
 636 const float Siox::FOREGROUND_CONFIDENCE=0.8f;
 637
 638 /**
 639  * Confidence for foreground or background type being equally likely.
 640  */
 641 const float Siox::UNKNOWN_REGION_CONFIDENCE=0.5f;
 642
 643 /**
 644  * Confidence for a region likely being background.
 645  */
 646 const float Siox::BACKGROUND_CONFIDENCE=0.1f;
 647
 648 /**
 649  * Confidence corresponding to a certain background reagion (equals zero).
 650  */
 651 const float Siox::CERTAIN_BACKGROUND_CONFIDENCE=0.0f;
 652
 653 /**
 654  *  Construct a Siox engine
 655  */
 656 Siox::Siox()
 657 {
 658     init();
 659 }
 660
 661
 662 /**
 663  *
 664  */
 665 Siox::~Siox()
 666 {
 667     cleanup();
 668 }
 669
 670
 671 /**
 672  * Error logging
 673  */
 674 void Siox::error(char *fmt, ...)
 675 {
 676     char msgbuf[256];
 677     va_list args;
 678     va_start(args, fmt);
 679     vsnprintf(msgbuf, 255, fmt, args);
 680     va_end(args) ;
 681 #ifdef HAVE_GLIB
 682     g_warning("Siox error: %s\n", msgbuf);
 683 #else
 684     fprintf(stderr, "Siox error: %s\n", msgbuf);
 685 #endif
 686 }
 687
 688 /**
 689  * Trace logging
 690  */
 691 void Siox::trace(char *fmt, ...)
 692 {
 693     char msgbuf[256];
 694     va_list args;
 695     va_start(args, fmt);
 696     vsnprintf(msgbuf, 255, fmt, args);
 697     va_end(args) ;
 698 #ifdef HAVE_GLIB
 699     g_message("Siox: %s\n", msgbuf);
 700 #else
 701     fprintf(stdout, "Siox: %s\n", msgbuf);
 702 #endif
 703 }
 704
 705
 706
 707
 708 /**
 709  *  Extract the foreground of the original image, according
 710  *  to the values in the confidence matrix.
 711  */
 712 SioxImage Siox::extractForeground(const SioxImage &originalImage,
 713                                   unsigned int backgroundFillColor)
 714 {
 715     init();
 716
 717     SioxImage workImage = originalImage;
 718
 719     //# fetch some info from the image
 720     width      = workImage.getWidth();
 721     height     = workImage.getHeight();
 722     pixelCount = width * height;
 723     image      = workImage.getImageData();
 724     cm         = workImage.getConfidenceData();
 725     labelField = new int[pixelCount];
 726
 727     trace("### Creating signatures");
 728
 729     //#### create color signatures
 730     std::vector<CLAB> knownBg;
 731     std::vector<CLAB> knownFg;
 732     for (int x = 0 ; x < workImage.getWidth() ; x++)
 733         for (int y = 0 ; y < workImage.getHeight() ; y++)
 734             {
 735             float cm = workImage.getConfidence(x, y);
 736             unsigned int pix = workImage.getPixel(x, y);
 737             if (cm <= BACKGROUND_CONFIDENCE)
 738                 knownBg.push_back(pix); //note: uses CLAB(rgb)
 739             else if (cm >= FOREGROUND_CONFIDENCE)
 740                 knownFg.push_back(pix);
 741             }
 742
 743     trace("knownBg:%d knownFg:%d", knownBg.size(), knownFg.size());
 744
 745
 746     std::vector<CLAB> bgSignature ;
 747     if (!colorSignature(knownBg, bgSignature, 3))
 748         {
 749         error("Could not create background signature");
 750         workImage.setValid(false);
 751         return workImage;
 752         }
 753     std::vector<CLAB> fgSignature ;
 754     if (!colorSignature(knownFg, fgSignature, 3))
 755         {
 756         error("Could not create foreground signature");
 757         delete[] labelField;
 758         workImage.setValid(false);
 759         return workImage;
 760         }
 761
 762     //trace("### bgSignature:%d", bgSignature.size());
 763
 764     if (bgSignature.size() < 1)
 765         {
 766         // segmentation impossible
 767         error("Signature size is < 1.  Segmentation is impossible");
 768         delete[] labelField;
 769         workImage.setValid(false);
 770         return workImage;
 771         }
 772
 773     // classify using color signatures,
 774     // classification cached in hashmap for drb and speedup purposes
 775     std::map<unsigned int, Tupel> hs;
 776
 777     for (unsigned int i=0; i<pixelCount; i++)
 778         {
 779         if (cm[i] >= FOREGROUND_CONFIDENCE)
 780             {
 781             cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
 782             }
 783         else if (cm[i] <= BACKGROUND_CONFIDENCE)
 784             {
 785             cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
 786             }
 787         else // somewhere in between
 788             {
 789             bool isBackground = true;
 790             std::map<unsigned int, Tupel>::iterator iter = hs.find(i);
 791             if (iter != hs.end()) //found
 792                 {
 793                 Tupel tupel = iter->second;
 794                 isBackground = tupel.minBgDist <= tupel.minFgDist;
 795                 }
 796             else
 797                 {
 798                 CLAB lab(image[i]);
 799                 float minBg = sqrEuclidianDist(lab, bgSignature[0]);
 800                 int minIndex=0;
 801                 for (unsigned int j=1; j<bgSignature.size() ; j++)
 802                     {
 803                     float d = sqrEuclidianDist(lab, bgSignature[j]);
 804                     if (d<minBg)
 805                         {
 806                         minBg    = d;
 807                         minIndex = j;
 808                         }
 809                     }
 810                 Tupel tupel(0.0f, 0,  0.0f, 0);
 811                 tupel.minBgDist  = minBg;
 812                 tupel.indexMinBg = minIndex;
 813                 float minFg = 1.0e6f;
 814                 minIndex = -1;
 815                 for (unsigned int j = 0 ; j < fgSignature.size() ; j++)
 816                     {
 817                     float d = sqrEuclidianDist(lab, fgSignature[j]);
 818                     if (d < minFg)
 819                         {
 820                         minFg    = d;
 821                         minIndex = j;
 822                         }
 823                     }
 824                 tupel.minFgDist  = minFg;
 825                 tupel.indexMinFg = minIndex;
 826                 if (fgSignature.size() == 0)
 827                     {
 828                     isBackground=(minBg <= clusterSize);
 829                     // remove next line to force behaviour of old algorithm
 830                     //error("foreground signature does not exist");
 831                     //delete[] labelField;
 832                     //workImage.setValid(false);
 833                     //return workImage;
 834                     }
 835                 else
 836                     {
 837                     isBackground = minBg < minFg;
 838                     }
 839                 hs[image[i]] = tupel;
 840                 }
 841
 842             if (isBackground)
 843                 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
 844             else
 845                 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
 846             }
 847         }
 848
 849     trace("### postProcessing");
 850
 851     //## postprocessing
 852     smooth(cm, width, height, 0.33f, 0.33f, 0.33f); // average
 853     normalizeMatrix(cm, pixelCount);
 854     erode(cm, width, height);
 855     keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.0/*sizeFactorToKeep*/);
 856
 857     for (int i=0; i < 2/*smoothness*/; i++)
 858         smooth(cm, width, height, 0.33f, 0.33f, 0.33f); // average
 859
 860     normalizeMatrix(cm, pixelCount);
 861
 862     for (unsigned int i=0; i < pixelCount; i++)
 863         {
 864         if (cm[i] >= UNKNOWN_REGION_CONFIDENCE)
 865             cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
 866         else
 867             cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
 868         }
 869
 870     keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.5/*sizeFactorToKeep*/);
 871     fillColorRegions();
 872     dilate(cm, width, height);
 873
 874     delete[] labelField;
 875
 876     //#### Yaay.  We are done.  Now clear everything but the background
 877     for (unsigned int y = 0 ; y < height ; y++)
 878         for (unsigned int x = 0 ; x < width ; x++)
 879             {
 880             float conf = workImage.getConfidence(x, y);
 881             if (conf < FOREGROUND_CONFIDENCE)
 882                 {
 883                 workImage.setPixel(x, y, backgroundFillColor);
 884                 }
 885             }
 886
 887     trace("### Done");
 888     return workImage;
 889 }
 890
 891
 892
 893 //##############
 894 //## PRIVATE
 895 //##############
 896
 897 /**
 898  *  Initialize the Siox engine to its 'pristine' state.
 899  *  Performed at the beginning of extractForeground().
 900  */
 901 void Siox::init()
 902 {
 903     limits[0] = 0.64f;
 904     limits[1] = 1.28f;
 905     limits[2] = 2.56f;
 906
 907     float negLimits[3];
 908     negLimits[0] = -limits[0];
 909     negLimits[1] = -limits[1];
 910     negLimits[2] = -limits[2];
 911
 912     clusterSize = sqrEuclidianDist(limits, 3, negLimits);
 913 }
 914
 915
 916 /**
 917  *  Clean up any debris from processing.
 918  */
 919 void Siox::cleanup()
 920 {
 921 }
 922
 923
 924
 925
 926 /**
 927  *  Stage 1 of the color signature work.  'dims' will be either
 928  *  2 for grays, or 3 for colors
 929  */
 930 void Siox::colorSignatureStage1(CLAB *points,
 931                                 unsigned int leftBase,
 932                                 unsigned int rightBase,
 933                                 unsigned int recursionDepth,
 934                                 unsigned int *clusterCount,
 935                                 const unsigned int dims)
 936 {
 937
 938     unsigned int currentDim = recursionDepth % dims;
 939     CLAB point = points[leftBase];
 940     float min = point(currentDim);
 941     float max = min;
 942
 943     for (unsigned int i = leftBase + 1; i < rightBase ; i++)
 944         {
 945         point = points[i];
 946         float curval = point(currentDim);
 947         if (curval < min) min = curval;
 948         if (curval > max) max = curval;
 949         }
 950
 951     //Do the Rubner-rule split (sounds like a dance)
 952     if (max - min > limits[currentDim])
 953         {
 954         float pivotPoint = (min + max) / 2.0; //average
 955         unsigned int left  = leftBase;
 956         unsigned int right = rightBase - 1;
 957
 958         //# partition points according to the dimension
 959         while (true)
 960             {
 961             while ( true )
 962                 {
 963                 point = points[left];
 964                 if (point(currentDim) > pivotPoint)
 965                     break;
 966                 left++;
 967                 }
 968             while ( true )
 969                 {
 970                 point = points[right];
 971                 if (point(currentDim) <= pivotPoint)
 972                     break;
 973                 right--;
 974                 }
 975
 976             if (left > right)
 977                 break;
 978
 979             point = points[left];
 980             points[left] = points[right];
 981             points[right] = point;
 982
 983             left++;
 984             right--;
 985             }
 986
 987         //# Recurse and create sub-trees
 988         colorSignatureStage1(points, leftBase, left,
 989                  recursionDepth + 1, clusterCount, dims);
 990         colorSignatureStage1(points, left, rightBase,
 991                  recursionDepth + 1, clusterCount, dims);
 992         }
 993     else
 994         {
 995         //create a leaf
 996         CLAB newpoint;
 997
 998         newpoint.C = rightBase - leftBase;
 999
1000         for (; leftBase < rightBase ; leftBase++)
1001             {
1002             newpoint.add(points[leftBase]);
1003             }
1004
1005         //printf("clusters:%d\n", *clusters);
1006
1007         if (newpoint.C != 0)
1008             newpoint.mul(1.0 / (float)newpoint.C);
1009         points[*clusterCount] = newpoint;
1010         (*clusterCount)++;
1011         }
1012 }
1013
1014
1015
1016 /**
1017  *  Stage 2 of the color signature work
1018  */
1019 void Siox::colorSignatureStage2(CLAB         *points,
1020                                 unsigned int leftBase,
1021                                 unsigned int rightBase,
1022                                 unsigned int recursionDepth,
1023                                 unsigned int *clusterCount,
1024                                 const float  threshold,
1025                                 const unsigned int dims)
1026 {
1027
1028
1029     unsigned int currentDim = recursionDepth % dims;
1030     CLAB point = points[leftBase];
1031     float min = point(currentDim);
1032     float max = min;
1033
1034     for (unsigned int i = leftBase+ 1; i < rightBase; i++)
1035         {
1036         point = points[i];
1037         float curval = point(currentDim);
1038         if (curval < min) min = curval;
1039         if (curval > max) max = curval;
1040         }
1041
1042     //Do the Rubner-rule split (sounds like a dance)
1043     if (max - min > limits[currentDim])
1044         {
1045         float pivotPoint = (min + max) / 2.0; //average
1046         unsigned int left  = leftBase;
1047         unsigned int right = rightBase - 1;
1048
1049         //# partition points according to the dimension
1050         while (true)
1051             {
1052             while ( true )
1053                 {
1054                 point = points[left];
1055                 if (point(currentDim) > pivotPoint)
1056                     break;
1057                 left++;
1058                 }
1059             while ( true )
1060                 {
1061                 point = points[right];
1062                 if (point(currentDim) <= pivotPoint)
1063                     break;
1064                 right--;
1065                 }
1066
1067             if (left > right)
1068                 break;
1069
1070             point = points[left];
1071             points[left] = points[right];
1072             points[right] = point;
1073
1074             left++;
1075             right--;
1076             }
1077
1078         //# Recurse and create sub-trees
1079         colorSignatureStage2(points, leftBase, left,
1080                  recursionDepth + 1, clusterCount, threshold, dims);
1081         colorSignatureStage2(points, left, rightBase,
1082                  recursionDepth + 1, clusterCount, threshold, dims);
1083         }
1084     else
1085         {
1086         //### Create a leaf
1087         unsigned int sum = 0;
1088         for (unsigned int i = leftBase; i < rightBase; i++)
1089             sum += points[i].C;
1090
1091         if ((float)sum >= threshold)
1092             {
1093             float scale = (float)(rightBase - leftBase);
1094             CLAB newpoint;
1095
1096             for (; leftBase < rightBase; leftBase++)
1097                 newpoint.add(points[leftBase]);
1098
1099             if (scale != 0.0)
1100                 newpoint.mul(1.0 / scale);
1101             points[*clusterCount] = newpoint;
1102             (*clusterCount)++;
1103             }
1104       }
1105 }
1106
1107
1108
1109 /**
1110  *  Main color signature method
1111  */
1112 bool Siox::colorSignature(const std::vector<CLAB> &inputVec,
1113                           std::vector<CLAB> &result,
1114                           const unsigned int dims)
1115 {
1116
1117     unsigned int length = inputVec.size();
1118
1119     if (length < 1) // no error. just don't do anything
1120         return true;
1121
1122     CLAB *input = (CLAB *) malloc(length * sizeof(CLAB));
1123
1124     if (!input)
1125         {
1126         error("Could not allocate buffer for signature");
1127         return false;
1128         }
1129     for (unsigned int i=0 ; i < length ; i++)
1130         input[i] = inputVec[i];
1131
1132     unsigned int stage1length = 0;
1133     colorSignatureStage1(input, 0, length, 0,
1134                    &stage1length, dims);
1135
1136     unsigned int stage2length = 0;
1137     colorSignatureStage2(input, 0, stage1length, 0,
1138                   &stage2length, length * 0.001, dims);
1139
1140     result.clear();
1141     for (unsigned int i=0 ; i < stage2length ; i++)
1142         result.push_back(input[i]);
1143
1144     free(input);
1145
1146     return true;
1147 }
1148
1149
1150
1151 /**
1152  *
1153  */
1154 void Siox::keepOnlyLargeComponents(float threshold,
1155                                    double sizeFactorToKeep)
1156 {
1157     for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1158         labelField[idx] = -1;
1159
1160     int curlabel = 0;
1161     int maxregion= 0;
1162     int maxblob  = 0;
1163
1164     // slow but easy to understand:
1165     std::vector<int> labelSizes;
1166     for (unsigned long i=0 ; i<pixelCount ; i++)
1167         {
1168         int regionCount = 0;
1169         if (labelField[i] == -1 && cm[i] >= threshold)
1170             {
1171             regionCount = depthFirstSearch(i, threshold, curlabel++);
1172             labelSizes.push_back(regionCount);
1173             }
1174
1175         if (regionCount>maxregion)
1176             {
1177             maxregion = regionCount;
1178             maxblob   = curlabel-1;
1179             }
1180         }
1181
1182     for (unsigned int i=0 ; i<pixelCount ; i++)
1183         {
1184         if (labelField[i] != -1)
1185             {
1186             // remove if the component is to small
1187             if (labelSizes[labelField[i]] * sizeFactorToKeep < maxregion)
1188                 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
1189
1190             // add maxblob always to foreground
1191             if (labelField[i] == maxblob)
1192                 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1193             }
1194         }
1195
1196 }
1197
1198
1199
1200 int Siox::depthFirstSearch(int startPos,
1201               float threshold, int curLabel)
1202 {
1203     // stores positions of labeled pixels, where the neighbours
1204     // should still be checked for processing:
1205
1206     //trace("startPos:%d threshold:%f curLabel:%d",
1207     //     startPos, threshold, curLabel);
1208
1209     std::vector<int> pixelsToVisit;
1210     int componentSize = 0;
1211
1212     if (labelField[startPos]==-1 && cm[startPos]>=threshold)
1213         {
1214         labelField[startPos] = curLabel;
1215         componentSize++;
1216         pixelsToVisit.push_back(startPos);
1217         }
1218
1219
1220     while (pixelsToVisit.size() > 0)
1221         {
1222         int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1223         pixelsToVisit.erase(pixelsToVisit.end() - 1);
1224         unsigned int x = pos % width;
1225         unsigned int y = pos / width;
1226
1227         // check all four neighbours
1228         int left = pos - 1;
1229         if (((int)x)-1>=0 && labelField[left]==-1 && cm[left]>=threshold)
1230             {
1231             labelField[left]=curLabel;
1232             componentSize++;
1233             pixelsToVisit.push_back(left);
1234             }
1235
1236         int right = pos + 1;
1237         if (x+1 < width && labelField[right]==-1 && cm[right]>=threshold)
1238             {
1239             labelField[right]=curLabel;
1240             componentSize++;
1241             pixelsToVisit.push_back(right);
1242             }
1243
1244         int top = pos - width;
1245         if (((int)y)-1>=0 && labelField[top]==-1 && cm[top]>=threshold)
1246             {
1247             labelField[top]=curLabel;
1248             componentSize++;
1249             pixelsToVisit.push_back(top);
1250             }
1251
1252         int bottom = pos + width;
1253         if (y+1 < height && labelField[bottom]==-1
1254                          && cm[bottom]>=threshold)
1255             {
1256             labelField[bottom]=curLabel;
1257             componentSize++;
1258             pixelsToVisit.push_back(bottom);
1259             }
1260
1261         }
1262     return componentSize;
1263 }
1264
1265
1266
1267 /**
1268  *
1269  */
1270 void Siox::fillColorRegions()
1271 {
1272     for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1273         labelField[idx] = -1;
1274
1275     //int maxRegion=0; // unused now
1276     std::vector<int> pixelsToVisit;
1277     for (unsigned long i=0; i<pixelCount; i++)
1278         { // for all pixels
1279         if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE)
1280             {
1281             continue; // already visited or bg
1282             }
1283
1284         unsigned int  origColor = image[i];
1285         unsigned long curLabel  = i+1;
1286         labelField[i]           = curLabel;
1287         cm[i]                   = CERTAIN_FOREGROUND_CONFIDENCE;
1288
1289         // int componentSize = 1;
1290         pixelsToVisit.push_back(i);
1291         // depth first search to fill region
1292         while (pixelsToVisit.size() > 0)
1293             {
1294             int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1295             pixelsToVisit.erase(pixelsToVisit.end() - 1);
1296             unsigned int x=pos % width;
1297             unsigned int y=pos / width;
1298             // check all four neighbours
1299             int left = pos-1;
1300             if (((int)x)-1 >= 0 && labelField[left] == -1
1301                         && CLAB::diff(image[left], origColor)<1.0)
1302                 {
1303                 labelField[left]=curLabel;
1304                 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1305                 // ++componentSize;
1306                 pixelsToVisit.push_back(left);
1307                 }
1308             int right = pos+1;
1309             if (x+1 < width && labelField[right]==-1
1310                         && CLAB::diff(image[right], origColor)<1.0)
1311                 {
1312                 labelField[right]=curLabel;
1313                 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1314                 // ++componentSize;
1315                 pixelsToVisit.push_back(right);
1316                 }
1317             int top = pos - width;
1318             if (((int)y)-1>=0 && labelField[top]==-1
1319                         && CLAB::diff(image[top], origColor)<1.0)
1320                 {
1321                 labelField[top]=curLabel;
1322                 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1323                 // ++componentSize;
1324                 pixelsToVisit.push_back(top);
1325                 }
1326             int bottom = pos + width;
1327             if (y+1 < height && labelField[bottom]==-1
1328                         && CLAB::diff(image[bottom], origColor)<1.0)
1329                 {
1330                 labelField[bottom]=curLabel;
1331                 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1332                 // ++componentSize;
1333                 pixelsToVisit.push_back(bottom);
1334                 }
1335             }
1336         //if (componentSize>maxRegion) {
1337         //    maxRegion=componentSize;
1338         //}
1339         }
1340 }
1341
1342
1343
1344
1345 /**
1346  * Applies the morphological dilate operator.
1347  *
1348  * Can be used to close small holes in the given confidence matrix.
1349  */
1350 void Siox::dilate(float *cm, int xres, int yres)
1351 {
1352
1353     for (int y=0; y<yres; y++)
1354         {
1355         for (int x=0; x<xres-1; x++)
1356              {
1357              int idx=(y*xres)+x;
1358              if (cm[idx+1]>cm[idx])
1359                  cm[idx]=cm[idx+1];
1360              }
1361         }
1362
1363     for (int y=0; y<yres; y++)
1364         {
1365         for (int x=xres-1; x>=1; x--)
1366             {
1367             int idx=(y*xres)+x;
1368             if (cm[idx-1]>cm[idx])
1369                 cm[idx]=cm[idx-1];
1370             }
1371         }
1372
1373     for (int y=0; y<yres-1; y++)
1374         {
1375         for (int x=0; x<xres; x++)
1376             {
1377             int idx=(y*xres)+x;
1378             if (cm[((y+1)*xres)+x] > cm[idx])
1379                 cm[idx]=cm[((y+1)*xres)+x];
1380             }
1381         }
1382
1383     for (int y=yres-1; y>=1; y--)
1384         {
1385         for (int x=0; x<xres; x++)
1386             {
1387             int idx=(y*xres)+x;
1388             if (cm[((y-1)*xres)+x] > cm[idx])
1389                 cm[idx]=cm[((y-1)*xres)+x];
1390             }
1391         }
1392 }
1393
1394 /**
1395  * Applies the morphological erode operator.
1396  */
1397 void Siox::erode(float *cm, int xres, int yres)
1398 {
1399     for (int y=0; y<yres; y++)
1400         {
1401         for (int x=0; x<xres-1; x++)
1402             {
1403             int idx=(y*xres)+x;
1404             if (cm[idx+1] < cm[idx])
1405                 cm[idx]=cm[idx+1];
1406             }
1407         }
1408     for (int y=0; y<yres; y++)
1409         {
1410         for (int x=xres-1; x>=1; x--)
1411             {
1412             int idx=(y*xres)+x;
1413             if (cm[idx-1] < cm[idx])
1414                 cm[idx]=cm[idx-1];
1415             }
1416         }
1417     for (int y=0; y<yres-1; y++)
1418         {
1419         for (int x=0; x<xres; x++)
1420             {
1421             int idx=(y*xres)+x;
1422             if (cm[((y+1)*xres)+x] < cm[idx])
1423                 cm[idx]=cm[((y+1)*xres)+x];
1424             }
1425         }
1426     for (int y=yres-1; y>=1; y--)
1427         {
1428         for (int x=0; x<xres; x++)
1429             {
1430             int idx=(y*xres)+x;
1431             if (cm[((y-1)*xres)+x] < cm[idx])
1432                 cm[idx]=cm[((y-1)*xres)+x];
1433             }
1434         }
1435 }
1436
1437
1438
1439
1440 /**
1441  * Normalizes the matrix to values to [0..1].
1442  */
1443 void Siox::normalizeMatrix(float *cm, int cmSize)
1444 {
1445     float max= -1000000.0f;
1446     for (int i=0; i<cmSize; i++)
1447         if (max<cm[i] > max)
1448             max=cm[i];
1449
1450     if (max<=0.0 || max==1.0)
1451         return;
1452
1453     float alpha=1.00f/max;
1454     premultiplyMatrix(alpha, cm, cmSize);
1455 }
1456
1457 /**
1458  * Multiplies matrix with the given scalar.
1459  */
1460 void Siox::premultiplyMatrix(float alpha, float *cm, int cmSize)
1461 {
1462     for (int i=0; i<cmSize; i++)
1463         cm[i]=alpha*cm[i];
1464 }
1465
1466 /**
1467  * Blurs confidence matrix with a given symmetrically weighted kernel.
1468  *
1469  * In the standard case confidence matrix entries are between 0...1 and
1470  * the weight factors sum up to 1.
1471  */
1472 void Siox::smooth(float *cm, int xres, int yres,
1473                   float f1, float f2, float f3)
1474 {
1475     for (int y=0; y<yres; y++)
1476         {
1477         for (int x=0; x<xres-2; x++)
1478             {
1479             int idx=(y*xres)+x;
1480             cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
1481             }
1482         }
1483     for (int y=0; y<yres; y++)
1484         {
1485         for (int x=xres-1; x>=2; x--)
1486             {
1487             int idx=(y*xres)+x;
1488             cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
1489             }
1490         }
1491     for (int y=0; y<yres-2; y++)
1492         {
1493         for (int x=0; x<xres; x++)
1494             {
1495             int idx=(y*xres)+x;
1496             cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
1497             }
1498         }
1499     for (int y=yres-1; y>=2; y--)
1500         {
1501         for (int x=0; x<xres; x++)
1502             {
1503             int idx=(y*xres)+x;
1504             cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
1505             }
1506         }
1507 }
1508
1509 /**
1510  * Squared Euclidian distance of p and q.
1511  */
1512 float Siox::sqrEuclidianDist(float *p, int pSize, float *q)
1513 {
1514     float sum=0.0;
1515     for (int i=0; i<pSize; i++)
1516         {
1517         float v = p[i] - q[i];
1518         sum += v*v;
1519         }
1520     return sum;
1521 }
1522
1523 /**
1524  * Squared Euclidian distance of p and q.
1525  */
1526 float Siox::sqrEuclidianDist(const CLAB &p, const CLAB &q)
1527 {
1528     float sum=0;
1529     sum += (p.L - q.L) * (p.L - q.L);
1530     sum += (p.A - q.A) * (p.A - q.A);
1531     sum += (p.B - q.B) * (p.B - q.B);
1532     return sum;
1533 }
1534
1535
1536
1537
1538
1539
1540
1541
1542 } // namespace siox
1543 } // namespace org
1544
1545 //########################################################################
1546 //#  E N D    O F    F I L E
1547 //########################################################################
1548