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(widthArg, heightArg);
 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(const 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(const 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(const 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     unsigned int currentDim = recursionDepth % dims;
1227     CieLab point = points[leftBase];
1228     float min = point(currentDim);
1229     float max = min;
1230
1231     for (unsigned int i = leftBase+ 1; i < rightBase; i++)
1232         {
1233         point = points[i];
1234         float curval = point(currentDim);
1235         if (curval < min) min = curval;
1236         if (curval > max) max = curval;
1237         }
1238
1239     //Do the Rubner-rule split (sounds like a dance)
1240     if (max - min > limits[currentDim])
1241         {
1242         float pivotPoint = (min + max) / 2.0; //average
1243         unsigned int left  = leftBase;
1244         unsigned int right = rightBase - 1;
1245
1246         //# partition points according to the dimension
1247         while (true)
1248             {
1249             while ( true )
1250                 {
1251                 point = points[left];
1252                 if (point(currentDim) > pivotPoint)
1253                     break;
1254                 left++;
1255                 }
1256             while ( true )
1257                 {
1258                 point = points[right];
1259                 if (point(currentDim) <= pivotPoint)
1260                     break;
1261                 right--;
1262                 }
1263
1264             if (left > right)
1265                 break;
1266
1267             point = points[left];
1268             points[left] = points[right];
1269             points[right] = point;
1270
1271             left++;
1272             right--;
1273             }
1274
1275         //# Recurse and create sub-trees
1276         colorSignatureStage2(points, leftBase, left,
1277                  recursionDepth + 1, clusterCount, threshold, dims);
1278         colorSignatureStage2(points, left, rightBase,
1279                  recursionDepth + 1, clusterCount, threshold, dims);
1280         }
1281     else
1282         {
1283         //### Create a leaf
1284         unsigned int sum = 0;
1285         for (unsigned int i = leftBase; i < rightBase; i++)
1286             sum += points[i].C;
1287
1288         if ((float)sum >= threshold)
1289             {
1290             float scale = (float)(rightBase - leftBase);
1291             CieLab newpoint;
1292
1293             for (; leftBase < rightBase; leftBase++)
1294                 newpoint.add(points[leftBase]);
1295
1296             if (scale != 0.0)
1297                 newpoint.mul(1.0 / scale);
1298             points[*clusterCount] = newpoint;
1299             (*clusterCount)++;
1300             }
1301       }
1302 }
1303
1304
1305
1306 /**
1307  *  Main color signature method
1308  */
1309 bool Siox::colorSignature(const std::vector<CieLab> &inputVec,
1310                           std::vector<CieLab> &result,
1311                           const unsigned int dims)
1312 {
1313
1314     unsigned int length = inputVec.size();
1315
1316     if (length < 1) // no error. just don't do anything
1317         return true;
1318
1319     CieLab *input = (CieLab *) malloc(length * sizeof(CieLab));
1320
1321     if (!input)
1322         {
1323         error("Could not allocate buffer for signature");
1324         return false;
1325         }
1326     for (unsigned int i=0 ; i < length ; i++)
1327         input[i] = inputVec[i];
1328
1329     unsigned int stage1length = 0;
1330     colorSignatureStage1(input, 0, length, 0,
1331                    &stage1length, dims);
1332
1333     unsigned int stage2length = 0;
1334     colorSignatureStage2(input, 0, stage1length, 0,
1335                   &stage2length, length * 0.001, dims);
1336
1337     result.clear();
1338     for (unsigned int i=0 ; i < stage2length ; i++)
1339         result.push_back(input[i]);
1340
1341     free(input);
1342
1343     return true;
1344 }
1345
1346
1347
1348 /**
1349  *
1350  */
1351 void Siox::keepOnlyLargeComponents(float threshold,
1352                                    double sizeFactorToKeep)
1353 {
1354     for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1355         labelField[idx] = -1;
1356
1357     int curlabel = 0;
1358     int maxregion= 0;
1359     int maxblob  = 0;
1360
1361     // slow but easy to understand:
1362     std::vector<int> labelSizes;
1363     for (unsigned long i=0 ; i<pixelCount ; i++)
1364         {
1365         int regionCount = 0;
1366         if (labelField[i] == -1 && cm[i] >= threshold)
1367             {
1368             regionCount = depthFirstSearch(i, threshold, curlabel++);
1369             labelSizes.push_back(regionCount);
1370             }
1371
1372         if (regionCount>maxregion)
1373             {
1374             maxregion = regionCount;
1375             maxblob   = curlabel-1;
1376             }
1377         }
1378
1379     for (unsigned int i=0 ; i<pixelCount ; i++)
1380         {
1381         if (labelField[i] != -1)
1382             {
1383             // remove if the component is to small
1384             if (labelSizes[labelField[i]] * sizeFactorToKeep < maxregion)
1385                 cm[i] = CERTAIN_BACKGROUND_CONFIDENCE;
1386
1387             // add maxblob always to foreground
1388             if (labelField[i] == maxblob)
1389                 cm[i] = CERTAIN_FOREGROUND_CONFIDENCE;
1390             }
1391         }
1392
1393 }
1394
1395
1396
1397 int Siox::depthFirstSearch(int startPos,
1398               float threshold, int curLabel)
1399 {
1400     // stores positions of labeled pixels, where the neighbours
1401     // should still be checked for processing:
1402
1403     //trace("startPos:%d threshold:%f curLabel:%d",
1404     //     startPos, threshold, curLabel);
1405
1406     std::vector<int> pixelsToVisit;
1407     int componentSize = 0;
1408
1409     if (labelField[startPos]==-1 && cm[startPos]>=threshold)
1410         {
1411         labelField[startPos] = curLabel;
1412         componentSize++;
1413         pixelsToVisit.push_back(startPos);
1414         }
1415
1416
1417     while (pixelsToVisit.size() > 0)
1418         {
1419         int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1420         pixelsToVisit.erase(pixelsToVisit.end() - 1);
1421         unsigned int x = pos % width;
1422         unsigned int y = pos / width;
1423
1424         // check all four neighbours
1425         int left = pos - 1;
1426         if (((int)x)-1>=0 && labelField[left]==-1 && cm[left]>=threshold)
1427             {
1428             labelField[left]=curLabel;
1429             componentSize++;
1430             pixelsToVisit.push_back(left);
1431             }
1432
1433         int right = pos + 1;
1434         if (x+1 < width && labelField[right]==-1 && cm[right]>=threshold)
1435             {
1436             labelField[right]=curLabel;
1437             componentSize++;
1438             pixelsToVisit.push_back(right);
1439             }
1440
1441         int top = pos - width;
1442         if (((int)y)-1>=0 && labelField[top]==-1 && cm[top]>=threshold)
1443             {
1444             labelField[top]=curLabel;
1445             componentSize++;
1446             pixelsToVisit.push_back(top);
1447             }
1448
1449         int bottom = pos + width;
1450         if (y+1 < height && labelField[bottom]==-1
1451                          && cm[bottom]>=threshold)
1452             {
1453             labelField[bottom]=curLabel;
1454             componentSize++;
1455             pixelsToVisit.push_back(bottom);
1456             }
1457
1458         }
1459     return componentSize;
1460 }
1461
1462
1463
1464 /**
1465  *
1466  */
1467 void Siox::fillColorRegions()
1468 {
1469     for (unsigned long idx = 0 ; idx<pixelCount ; idx++)
1470         labelField[idx] = -1;
1471
1472     //int maxRegion=0; // unused now
1473     std::vector<int> pixelsToVisit;
1474     for (unsigned long i=0; i<pixelCount; i++)
1475         { // for all pixels
1476         if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE)
1477             {
1478             continue; // already visited or bg
1479             }
1480
1481         unsigned int  origColor = image[i];
1482         unsigned long curLabel  = i+1;
1483         labelField[i]           = curLabel;
1484         cm[i]                   = CERTAIN_FOREGROUND_CONFIDENCE;
1485
1486         // int componentSize = 1;
1487         pixelsToVisit.push_back(i);
1488         // depth first search to fill region
1489         while (pixelsToVisit.size() > 0)
1490             {
1491             int pos = pixelsToVisit[pixelsToVisit.size() - 1];
1492             pixelsToVisit.erase(pixelsToVisit.end() - 1);
1493             unsigned int x=pos % width;
1494             unsigned int y=pos / width;
1495             // check all four neighbours
1496             int left = pos-1;
1497             if (((int)x)-1 >= 0 && labelField[left] == -1
1498                         && CieLab::diff(image[left], origColor)<1.0)
1499                 {
1500                 labelField[left]=curLabel;
1501                 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1502                 // ++componentSize;
1503                 pixelsToVisit.push_back(left);
1504                 }
1505             int right = pos+1;
1506             if (x+1 < width && labelField[right]==-1
1507                         && CieLab::diff(image[right], origColor)<1.0)
1508                 {
1509                 labelField[right]=curLabel;
1510                 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1511                 // ++componentSize;
1512                 pixelsToVisit.push_back(right);
1513                 }
1514             int top = pos - width;
1515             if (((int)y)-1>=0 && labelField[top]==-1
1516                         && CieLab::diff(image[top], origColor)<1.0)
1517                 {
1518                 labelField[top]=curLabel;
1519                 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1520                 // ++componentSize;
1521                 pixelsToVisit.push_back(top);
1522                 }
1523             int bottom = pos + width;
1524             if (y+1 < height && labelField[bottom]==-1
1525                         && CieLab::diff(image[bottom], origColor)<1.0)
1526                 {
1527                 labelField[bottom]=curLabel;
1528                 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1529                 // ++componentSize;
1530                 pixelsToVisit.push_back(bottom);
1531                 }
1532             }
1533         //if (componentSize>maxRegion) {
1534         //    maxRegion=componentSize;
1535         //}
1536         }
1537 }
1538
1539
1540
1541
1542 /**
1543  * Applies the morphological dilate operator.
1544  *
1545  * Can be used to close small holes in the given confidence matrix.
1546  */
1547 void Siox::dilate(float *cm, int xres, int yres)
1548 {
1549
1550     for (int y=0; y<yres; y++)
1551         {
1552         for (int x=0; x<xres-1; x++)
1553              {
1554              int idx=(y*xres)+x;
1555              if (cm[idx+1]>cm[idx])
1556                  cm[idx]=cm[idx+1];
1557              }
1558         }
1559
1560     for (int y=0; y<yres; y++)
1561         {
1562         for (int x=xres-1; x>=1; x--)
1563             {
1564             int idx=(y*xres)+x;
1565             if (cm[idx-1]>cm[idx])
1566                 cm[idx]=cm[idx-1];
1567             }
1568         }
1569
1570     for (int y=0; y<yres-1; y++)
1571         {
1572         for (int x=0; x<xres; x++)
1573             {
1574             int idx=(y*xres)+x;
1575             if (cm[((y+1)*xres)+x] > cm[idx])
1576                 cm[idx]=cm[((y+1)*xres)+x];
1577             }
1578         }
1579
1580     for (int y=yres-1; y>=1; y--)
1581         {
1582         for (int x=0; x<xres; x++)
1583             {
1584             int idx=(y*xres)+x;
1585             if (cm[((y-1)*xres)+x] > cm[idx])
1586                 cm[idx]=cm[((y-1)*xres)+x];
1587             }
1588         }
1589 }
1590
1591 /**
1592  * Applies the morphological erode operator.
1593  */
1594 void Siox::erode(float *cm, int xres, int yres)
1595 {
1596     for (int y=0; y<yres; y++)
1597         {
1598         for (int x=0; x<xres-1; x++)
1599             {
1600             int idx=(y*xres)+x;
1601             if (cm[idx+1] < cm[idx])
1602                 cm[idx]=cm[idx+1];
1603             }
1604         }
1605     for (int y=0; y<yres; y++)
1606         {
1607         for (int x=xres-1; x>=1; x--)
1608             {
1609             int idx=(y*xres)+x;
1610             if (cm[idx-1] < cm[idx])
1611                 cm[idx]=cm[idx-1];
1612             }
1613         }
1614     for (int y=0; y<yres-1; y++)
1615         {
1616         for (int x=0; x<xres; x++)
1617             {
1618             int idx=(y*xres)+x;
1619             if (cm[((y+1)*xres)+x] < cm[idx])
1620                 cm[idx]=cm[((y+1)*xres)+x];
1621             }
1622         }
1623     for (int y=yres-1; y>=1; y--)
1624         {
1625         for (int x=0; x<xres; x++)
1626             {
1627             int idx=(y*xres)+x;
1628             if (cm[((y-1)*xres)+x] < cm[idx])
1629                 cm[idx]=cm[((y-1)*xres)+x];
1630             }
1631         }
1632 }
1633
1634
1635
1636 /**
1637  * Normalizes the matrix to values to [0..1].
1638  */
1639 void Siox::normalizeMatrix(float *cm, int cmSize)
1640 {
1641     float max= -1000000.0f;
1642     for (int i=0; i<cmSize; i++)
1643         if (cm[i] > max) max=cm[i];
1644
1645     //good to use STL, but max() is not iterative
1646     //float max = *std::max(cm, cm + cmSize);
1647
1648     if (max<=0.0 || max==1.0)
1649         return;
1650
1651     float alpha=1.00f/max;
1652     premultiplyMatrix(alpha, cm, cmSize);
1653 }
1654
1655 /**
1656  * Multiplies matrix with the given scalar.
1657  */
1658 void Siox::premultiplyMatrix(float alpha, float *cm, int cmSize)
1659 {
1660     for (int i=0; i<cmSize; i++)
1661         cm[i]=alpha*cm[i];
1662 }
1663
1664 /**
1665  * Blurs confidence matrix with a given symmetrically weighted kernel.
1666  *
1667  * In the standard case confidence matrix entries are between 0...1 and
1668  * the weight factors sum up to 1.
1669  */
1670 void Siox::smooth(float *cm, int xres, int yres,
1671                   float f1, float f2, float f3)
1672 {
1673     for (int y=0; y<yres; y++)
1674         {
1675         for (int x=0; x<xres-2; x++)
1676             {
1677             int idx=(y*xres)+x;
1678             cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
1679             }
1680         }
1681     for (int y=0; y<yres; y++)
1682         {
1683         for (int x=xres-1; x>=2; x--)
1684             {
1685             int idx=(y*xres)+x;
1686             cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
1687             }
1688         }
1689     for (int y=0; y<yres-2; y++)
1690         {
1691         for (int x=0; x<xres; x++)
1692             {
1693             int idx=(y*xres)+x;
1694             cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
1695             }
1696         }
1697     for (int y=yres-1; y>=2; y--)
1698         {
1699         for (int x=0; x<xres; x++)
1700             {
1701             int idx=(y*xres)+x;
1702             cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
1703             }
1704         }
1705 }
1706
1707 /**
1708  * Squared Euclidian distance of p and q.
1709  */
1710 float Siox::sqrEuclidianDist(float *p, int pSize, float *q)
1711 {
1712     float sum=0.0;
1713     for (int i=0; i<pSize; i++)
1714         {
1715         float v = p[i] - q[i];
1716         sum += v*v;
1717         }
1718     return sum;
1719 }
1720
1721
1722
1723
1724
1725
1726 } // namespace siox
1727 } // namespace org
1728
1729 //########################################################################
1730 //#  E N D    O F    F I L E
1731 //########################################################################
1732