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