src/trace/siox.cpp

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