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 <string>
  21
  22 #include <stdarg.h> //for error() and trace()
  23 #include <math.h>   //sqrt(), pow(), round(), etc
  24
  25
  26 namespace org
  27 {
  28 namespace siox
  29 {
  30
  31 //########################################################################
  32 //##  U T I L S    (originally Utils.java)
  33 //########################################################################
  34
  35 /**
  36  * Collection of auxiliary image processing methods used by the
  37  * SioxSegmentator mainly for postprocessing.
  38  *
  39  * @author G. Friedland, K. Jantz, L. Knipping
  40  * @version 1.05
  41  *
  42  * Conversion to C++ by Bob Jamison
  43  *
  44  */
  45
  46
  47 /** Caches color conversion values to speed up RGB->CIELAB conversion.*/
  48 static std::map<unsigned long, CLAB> RGB_TO_LAB;
  49
  50 //forward decls
  51 static void premultiplyMatrix(float alpha, float *cm, int cmSize);
  52 //static float colordiffsq(unsigned long rgb0, unsigned long rgb1);
  53 //static int getAlpha(unsigned long argb);
  54 static int getRed(unsigned long rgb);
  55 static int getGreen(unsigned long rgb);
  56 static int getBlue(unsigned long rgb);
  57 //static unsigned long packPixel(int a, int r, int g, int b);
  58 static CLAB rgbToClab(unsigned long rgb);
  59
  60 /**
  61  * Applies the morphological dilate operator.
  62  *
  63  * Can be used to close small holes in the given confidence matrix.
  64  *
  65  * @param cm Confidence matrix to be processed.
  66  * @param xres Horizontal resolution of the matrix.
  67  * @param yres Vertical resolution of the matrix.
  68  */
  69 static void dilate(float *cm, int xres, int yres)
  70 {
  71     for (int y=0; y<yres; y++) {
  72         for (int x=0; x<xres-1; x++) {
  73             int idx=(y*xres)+x;
  74             if (cm[idx+1]>cm[idx])
  75                 cm[idx]=cm[idx+1];
  76             }
  77         }
  78     for (int y=0; y<yres; y++) {
  79             for (int x=xres-1; x>=1; x--) {
  80             int idx=(y*xres)+x;
  81             if (cm[idx-1]>cm[idx])
  82                 cm[idx]=cm[idx-1];
  83         }
  84     }
  85     for (int y=0; y<yres-1; y++) {
  86         for (int x=0; x<xres; x++) {
  87         int idx=(y*xres)+x;
  88         if (cm[((y+1)*xres)+x] > cm[idx])
  89             cm[idx]=cm[((y+1)*xres)+x];
  90         }
  91     }
  92     for (int y=yres-1; y>=1; y--) {
  93         for (int x=0; x<xres; x++) {
  94         int idx=(y*xres)+x;
  95         if (cm[((y-1)*xres)+x] > cm[idx])
  96             cm[idx]=cm[((y-1)*xres)+x];
  97         }
  98     }
  99 }
 100
 101 /**
 102  * Applies the morphological erode operator.
 103  *
 104  * @param cm Confidence matrix to be processed.
 105  * @param xres Horizontal resolution of the matrix.
 106  * @param yres Vertical resolution of the matrix.
 107  */
 108 static void erode(float *cm, int xres, int yres)
 109 {
 110     for (int y=0; y<yres; y++) {
 111         for (int x=0; x<xres-1; x++) {
 112             int idx=(y*xres)+x;
 113             if (cm[idx+1] < cm[idx])
 114                 cm[idx]=cm[idx+1];
 115         }
 116     }
 117     for (int y=0; y<yres; y++) {
 118         for (int x=xres-1; x>=1; x--) {
 119             int idx=(y*xres)+x;
 120             if (cm[idx-1] < cm[idx])
 121                 cm[idx]=cm[idx-1];
 122         }
 123     }
 124     for (int y=0; y<yres-1; y++) {
 125         for (int x=0; x<xres; x++) {
 126             int idx=(y*xres)+x;
 127             if (cm[((y+1)*xres)+x] < cm[idx])
 128                 cm[idx]=cm[((y+1)*xres)+x];
 129         }
 130     }
 131     for (int y=yres-1; y>=1; y--) {
 132         for (int x=0; x<xres; x++) {
 133             int idx=(y*xres)+x;
 134             if (cm[((y-1)*xres)+x] < cm[idx])
 135                 cm[idx]=cm[((y-1)*xres)+x];
 136         }
 137     }
 138 }
 139
 140 /**
 141  * Normalizes the matrix to values to [0..1].
 142  *
 143  * @param cm The matrix to be normalized.
 144  */
 145 static void normalizeMatrix(float *cm, int cmSize)
 146 {
 147     float max=0.0f;
 148     for (int i=0; i<cmSize; i++) {
 149         if (max<cm[i])
 150             max=cm[i];
 151     }
 152     if (max<=0.0)
 153         return;
 154     else if (max==1.00)
 155         return;
 156
 157     float alpha=1.00f/max;
 158     premultiplyMatrix(alpha, cm, cmSize);
 159 }
 160
 161 /**
 162  * Multiplies matrix with the given scalar.
 163  *
 164  * @param alpha The scalar value.
 165  * @param cm The matrix of values be multiplied with alpha.
 166  * @param cmSize The matrix length.
 167  */
 168 static void premultiplyMatrix(float alpha, float *cm, int cmSize)
 169 {
 170     for (int i=0; i<cmSize; i++)
 171         cm[i]=alpha*cm[i];
 172 }
 173
 174 /**
 175  * Blurs confidence matrix with a given symmetrically weighted kernel.
 176  * <P>
 177  * In the standard case confidence matrix entries are between 0...1 and
 178  * the weight factors sum up to 1.
 179  *
 180  * @param cm The matrix to be smoothed.
 181  * @param xres Horizontal resolution of the matrix.
 182  * @param yres Vertical resolution of the matrix.
 183  * @param f1 Weight factor for the first pixel.
 184  * @param f2 Weight factor for the mid-pixel.
 185  * @param f3 Weight factor for the last pixel.
 186  */
 187 static void smoothcm(float *cm, int xres, int yres,
 188                      float f1, float f2, float f3)
 189 {
 190     for (int y=0; y<yres; y++) {
 191         for (int x=0; x<xres-2; x++) {
 192             int idx=(y*xres)+x;
 193             cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
 194         }
 195     }
 196     for (int y=0; y<yres; y++) {
 197         for (int x=xres-1; x>=2; x--) {
 198             int idx=(y*xres)+x;
 199             cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
 200         }
 201     }
 202     for (int y=0; y<yres-2; y++) {
 203         for (int x=0; x<xres; x++) {
 204             int idx=(y*xres)+x;
 205             cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
 206         }
 207     }
 208     for (int y=yres-1; y>=2; y--) {
 209         for (int x=0; x<xres; x++) {
 210             int idx=(y*xres)+x;
 211             cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
 212         }
 213     }
 214 }
 215
 216 /**
 217  * Squared Euclidian distance of p and q.
 218  * <P>
 219  * Usage hint: When only comparisons between Euclidian distances without
 220  * actual values are needed, the squared distance can be preferred
 221  * for being faster to compute.
 222  *
 223  * @param p First euclidian point coordinates.
 224  * @param pSize Length of coordinate array.
 225  * @param q Second euclidian point coordinates.
 226  *        Dimension must not be smaller than that of p.
 227  *        Any extra dimensions will be ignored.
 228  * @return Squared euclidian distance of p and q.
 229  * @see #euclid
 230  */
 231 static float sqrEuclidianDist(float *p, int pSize, float *q)
 232 {
 233     float sum=0;
 234     for (int i=0; i<pSize; i++)
 235         sum+=(p[i]-q[i])*(p[i]-q[i]);
 236     return sum;
 237 }
 238
 239 /**
 240  * Squared Euclidian distance of p and q.
 241  * <P>
 242  * Usage hint: When only comparisons between Euclidian distances without
 243  * actual values are needed, the squared distance can be preferred
 244  * for being faster to compute.
 245  *
 246  * @param p CLAB value
 247  * @param q second CLAB value
 248  * @return Squared euclidian distance of p and q.
 249  * @see #euclid
 250  */
 251 static float sqrEuclidianDist(const CLAB &p, const CLAB &q)
 252 {
 253     float sum=0;
 254     sum += (p.C - q.C) * (p.C - q.C);
 255     sum += (p.L - q.L) * (p.L - q.L);
 256     sum += (p.A - q.A) * (p.A - q.A);
 257     sum += (p.B - q.B) * (p.B - q.B);
 258     return sum;
 259 }
 260
 261 /**
 262  * Euclidian distance of p and q.
 263  *
 264  * @param p First euclidian point coordinates.
 265  * @param pSize Length of coordinate array.
 266  * @param q Second euclidian point coordinates.
 267  *        Dimension must not be smaller than that of p.
 268  *        Any extra dimensions will be ignored.
 269  * @return Squared euclidian distance of p and q.
 270  * @see #sqrEuclidianDist
 271  */
 272 /*
 273 static float euclid(float *p, int pSize, float *q)
 274 {
 275     return (float)sqrt(sqrEuclidianDist(p, pSize, q));
 276 }
 277 */
 278
 279 /**
 280  * Computes Euclidian distance of two RGB color values.
 281  *
 282  * @param rgb0 First color value.
 283  * @param rgb1 Second color value.
 284  * @return Euclidian distance between the two color values.
 285  */
 286 /*
 287 static float colordiff(long rgb0, long rgb1)
 288 {
 289     return (float)sqrt(colordiffsq(rgb0, rgb1));
 290 }
 291 */
 292
 293 /**
 294  * Computes squared euclidian distance of two RGB color values.
 295  * <P>
 296  * Note: Faster to compute than colordiff
 297  *
 298  * @param rgb0 First color value.
 299  * @param rgb1 Second color value.
 300  * @return Squared Euclidian distance between the two color values.
 301  */
 302 /*
 303 static float colordiffsq(long rgb0, long rgb1)
 304 {
 305     int rDist=getRed(rgb0)   - getRed(rgb1);
 306     int gDist=getGreen(rgb0) - getGreen(rgb1);
 307     int bDist=getBlue(rgb0)  - getBlue(rgb1);
 308
 309     return (float)(rDist*rDist+gDist*gDist+bDist*bDist);
 310 }
 311 */
 312
 313 /**
 314  * Averages two ARGB colors.
 315  *
 316  * @param argb0 First color value.
 317  * @param argb1 Second color value.
 318  * @return The averaged ARGB color.
 319  */
 320 /*
 321 static long average(long argb0, long argb1)
 322 {
 323     long ret = packPixel(
 324          (getAlpha(argb0) + getAlpha(argb1))/2,
 325          (getRed(argb0)   + getRed(argb1)  )/2,
 326          (getGreen(argb0) + getGreen(argb1))/2,
 327          (getBlue(argb0)  + getBlue(argb1) )/2);
 328     return ret;
 329 }
 330 */
 331
 332 /**
 333  * Computes squared euclidian distance in CLAB space for two colors
 334  * given as RGB values.
 335  *
 336  * @param rgb0 First color value.
 337  * @param rgb1 Second color value.
 338  * @return Squared Euclidian distance in CLAB space.
 339  */
 340 static float labcolordiffsq(unsigned long rgb1, unsigned long rgb2)
 341 {
 342     CLAB c1 = rgbToClab(rgb1);
 343     CLAB c2 = rgbToClab(rgb2);
 344     float euclid=0.0f;
 345     euclid += (c1.L - c2.L) * (c1.L - c2.L);
 346     euclid += (c1.A - c2.A) * (c1.A - c2.A);
 347     euclid += (c1.B - c2.B) * (c1.B - c2.B);
 348     return euclid;
 349 }
 350
 351
 352 /**
 353  * Computes squared euclidian distance in CLAB space for two colors
 354  * given as RGB values.
 355  *
 356  * @param rgb0 First color value.
 357  * @param rgb1 Second color value.
 358  * @return Euclidian distance in CLAB space.
 359  */
 360 static float labcolordiff(unsigned long rgb0, unsigned long rgb1)
 361 {
 362     return (float)sqrt(labcolordiffsq(rgb0, rgb1));
 363 }
 364
 365
 366 /**
 367  * Converts 24-bit RGB values to {l,a,b} float values.
 368  * <P>
 369  * The conversion used is decribed at
 370  * <a href="http://www.easyrgb.com/math.php?MATH=M7#text7">CLAB Conversion</a>
 371  * for reference white D65. Note that that the conversion is computational
 372  * expensive. Result are cached to speed up further conversion calls.
 373  *
 374  * @param rgb RGB color value,
 375  * @return CLAB color value tripel.
 376  */
 377 static CLAB rgbToClab(unsigned long rgb)
 378 {
 379     std::map<unsigned long, CLAB>::iterator iter = RGB_TO_LAB.find(rgb);
 380     if (iter != RGB_TO_LAB.end())
 381         {
 382         CLAB res = iter->second;
 383         return res;
 384         }
 385
 386     int R=getRed(rgb);
 387     int G=getGreen(rgb);
 388     int B=getBlue(rgb);
 389
 390     float var_R=(R/255.0f); //R = From 0 to 255
 391     float var_G=(G/255.0f); //G = From 0 to 255
 392     float var_B=(B/255.0f); //B = From 0 to 255
 393
 394     if (var_R>0.04045)
 395         var_R=(float) pow((var_R+0.055f)/1.055f, 2.4);
 396     else
 397         var_R=var_R/12.92f;
 398
 399     if (var_G>0.04045)
 400         var_G=(float) pow((var_G+0.055f)/1.055f, 2.4);
 401     else
 402         var_G=var_G/12.92f;
 403
 404     if (var_B>0.04045)
 405         var_B=(float) pow((var_B+0.055f)/1.055f, 2.4);
 406     else
 407         var_B=var_B/12.92f;
 408
 409     var_R=var_R*100.0f;
 410     var_G=var_G*100.0f;
 411     var_B=var_B*100.0f;
 412
 413     // Observer. = 2�, Illuminant = D65
 414     float X=var_R*0.4124f + var_G*0.3576f + var_B*0.1805f;
 415     float Y=var_R*0.2126f + var_G*0.7152f + var_B*0.0722f;
 416     float Z=var_R*0.0193f + var_G*0.1192f + var_B*0.9505f;
 417
 418     float var_X=X/95.047f;
 419     float var_Y=Y/100.0f;
 420     float var_Z=Z/108.883f;
 421
 422     if (var_X>0.008856f)
 423         var_X=(float) pow(var_X, 0.3333f);
 424     else
 425         var_X=(7.787f*var_X)+(16.0f/116.0f);
 426
 427     if (var_Y>0.008856f)
 428         var_Y=(float) pow(var_Y, 0.3333f);
 429     else
 430         var_Y=(7.787f*var_Y)+(16.0f/116.0f);
 431
 432     if (var_Z>0.008856f)
 433         var_Z=(float) pow(var_Z, 0.3333f);
 434     else
 435         var_Z=(7.787f*var_Z)+(16.0f/116.0f);
 436
 437     CLAB lab((116.0f*var_Y)-16.0f , 500.0f*(var_X-var_Y), 200.0f*(var_Y-var_Z));
 438
 439     RGB_TO_LAB[rgb] = lab;
 440
 441     return lab;
 442 }
 443
 444 /**
 445  * Converts an CLAB value to a RGB color value.
 446  * <P>
 447  * This is the reverse operation to rgbToClab.
 448  * @param clab CLAB value.
 449  * @return RGB value.
 450  * @see #rgbToClab
 451  */
 452 /*
 453 static long clabToRGB(const CLAB &clab)
 454 {
 455     float L=clab.L;
 456     float a=clab.A;
 457     float b=clab.B;
 458
 459     float var_Y=(L+16.0f)/116.0f;
 460     float var_X=a/500.0f+var_Y;
 461     float var_Z=var_Y-b/200.0f;
 462
 463     float var_yPow3=(float)pow(var_Y, 3.0);
 464     float var_xPow3=(float)pow(var_X, 3.0);
 465     float var_zPow3=(float)pow(var_Z, 3.0);
 466
 467     if (var_yPow3>0.008856f)
 468         var_Y=var_yPow3;
 469     else
 470         var_Y=(var_Y-16.0f/116.0f)/7.787f;
 471
 472     if (var_xPow3>0.008856f)
 473         var_X=var_xPow3;
 474     else
 475         var_X=(var_X-16.0f/116.0f)/7.787f;
 476
 477     if (var_zPow3>0.008856f)
 478         var_Z=var_zPow3;
 479     else
 480         var_Z=(var_Z-16.0f/116.0f)/7.787f;
 481
 482     float X=95.047f  * var_X; //ref_X= 95.047  Observer=2�, Illuminant=D65
 483     float Y=100.0f   * var_Y; //ref_Y=100.000
 484     float Z=108.883f * var_Z; //ref_Z=108.883
 485
 486     var_X=X/100.0f; //X = From 0 to ref_X
 487     var_Y=Y/100.0f; //Y = From 0 to ref_Y
 488     var_Z=Z/100.0f; //Z = From 0 to ref_Y
 489
 490     float var_R=(float)(var_X *  3.2406f + var_Y * -1.5372f + var_Z * -0.4986f);
 491     float var_G=(float)(var_X * -0.9689f + var_Y *  1.8758f + var_Z *  0.0415f);
 492     float var_B=(float)(var_X *  0.0557f + var_Y * -0.2040f + var_Z *  1.0570f);
 493
 494     if (var_R>0.0031308f)
 495         var_R=(float)(1.055f*pow(var_R, (1.0f/2.4f))-0.055f);
 496     else
 497         var_R=12.92f*var_R;
 498
 499     if (var_G>0.0031308f)
 500         var_G=(float)(1.055f*pow(var_G, (1.0f/2.4f))-0.055f);
 501     else
 502         var_G=12.92f*var_G;
 503
 504     if (var_B>0.0031308f)
 505         var_B=(float)(1.055f*pow(var_B, (1.0f/2.4f))-0.055f);
 506     else
 507         var_B=12.92f*var_B;
 508
 509     int R = (int)lround(var_R*255.0f);
 510     int G = (int)lround(var_G*255.0f);
 511     int B = (int)lround(var_B*255.0f);
 512
 513     return packPixel(0xFF, R, G, B);
 514 }
 515 */
 516
 517 /**
 518  * Sets the alpha byte of a pixel.
 519  *
 520  * Constructs alpha to values from 0 to 255.
 521  * @param alpha Alpha value from 0 (transparent) to 255 (opaque).
 522  * @param rgb The 24bit rgb color to be combined with the alpga value.
 523  * @return An ARBG calor value.
 524  */
 525 static long setAlpha(int alpha, unsigned long rgb)
 526 {
 527     if (alpha>255)
 528         alpha=0;
 529     else if (alpha<0)
 530         alpha=0;
 531     return (alpha<<24)|(rgb&0xFFFFFF);
 532 }
 533
 534 /**
 535  * Sets the alpha byte of a pixel.
 536  *
 537  * Constricts alpha to values from 0 to 255.
 538  * @param alpha Alpha value from 0.0f (transparent) to 1.0f (opaque).
 539  * @param rgb The 24bit rgb color to be combined with the alpga value.
 540  * @return An ARBG calor value.
 541  */
 542 static long setAlpha(float alpha, unsigned long rgb)
 543 {
 544     return setAlpha((int)(255.0f*alpha), rgb);
 545 }
 546
 547 /**
 548  * Limits the values of a,r,g,b to values from 0 to 255 and puts them
 549  * together into an 32 bit integer.
 550  *
 551  * @param a Alpha part, the first byte.
 552  * @param r Red part, the second byte.
 553  * @param g Green part, the third byte.
 554  * @param b Blue part, the fourth byte.
 555  * @return A ARBG value packed to an int.
 556  */
 557 /*
 558 static long packPixel(int a, int r, int g, int b)
 559 {
 560     if (a<0)
 561         a=0;
 562     else if (a>255)
 563         a=255;
 564
 565     if (r<0)
 566         r=0;
 567     else if (r>255)
 568         r=255;
 569
 570     if (g<0)
 571         g=0;
 572     else if (g>255)
 573         g=255;
 574
 575     if (b<0)
 576         b=0;
 577     else if (b>255)
 578         b=255;
 579
 580     return (a<<24)|(r<<16)|(g<<8)|b;
 581 }
 582 */
 583
 584 /**
 585  * Returns the alpha component of an ARGB value.
 586  *
 587  * @param argb An ARGB color value.
 588  * @return The alpha component, ranging from 0 to 255.
 589  */
 590 /*
 591 static int getAlpha(unsigned long argb)
 592 {
 593     return (argb>>24)&0xFF;
 594 }
 595 */
 596
 597 /**
 598  * Returns the red component of an (A)RGB value.
 599  *
 600  * @param rgb An (A)RGB color value.
 601  * @return The red component, ranging from 0 to 255.
 602  */
 603 static int getRed(unsigned long rgb)
 604 {
 605     return (rgb>>16)&0xFF;
 606 }
 607
 608
 609 /**
 610  * Returns the green component of an (A)RGB value.
 611  *
 612  * @param rgb An (A)RGB color value.
 613  * @return The green component, ranging from 0 to 255.
 614  */
 615 static int getGreen(unsigned long rgb)
 616 {
 617     return (rgb>>8)&0xFF;
 618 }
 619
 620 /**
 621  * Returns the blue component of an (A)RGB value.
 622  *
 623  * @param rgb An (A)RGB color value.
 624  * @return The blue component, ranging from 0 to 255.
 625  */
 626 static int getBlue(unsigned long rgb)
 627 {
 628     return (rgb)&0xFF;
 629 }
 630
 631 /**
 632  * Returns a string representation of a CLAB value.
 633  *
 634  * @param clab The CLAB value.
 635  * @param clabSize Size of the CLAB value.
 636  * @return A string representation of the CLAB value.
 637  */
 638 /*
 639 static std::string clabToString(const CLAB &clab)
 640 {
 641     std::string buff;
 642     char nbuf[60];
 643     snprintf(nbuf, 59, "%5.3f, %5.3f, %5.3f", clab.L, clab.A, clab.B);
 644     buff = nbuf;
 645     return buff;
 646 }
 647 */
 648
 649 //########################################################################
 650 //##  C O L O R    S I G N A T U R E    (originally ColorSignature.java)
 651 //########################################################################
 652
 653 /**
 654  * Representation of a color signature.
 655  * <br><br>
 656  * This class implements a clustering algorithm based on a modified kd-tree.
 657  * The splitting rule is to simply divide the given interval into two equally
 658  * sized subintervals.
 659  * In the <code>stageone()</code>, approximate clusters are found by building
 660  * up such a tree and stopping when an interval at a node has become smaller
 661  * than the allowed cluster diameter, which is given by <code>limits</code>.
 662  * At this point, clusters may be split in several nodes.<br>
 663  * Therefore, in <code>stagetwo()</code>, nodes that belong to several clusters
 664  * are recombined by another k-d tree clustering on the prior cluster
 665  * centroids. To guarantee a proper level of abstraction, clusters that contain
 666  * less than 0.01% of the pixels of the entire sample are removed. Please
 667  * refer to the file NOTICE to get links to further documentation.
 668  *
 669  * @author Gerald Friedland, Lars Knipping
 670  * @version 1.02
 671  *
 672  * Conversion to C++ by Bob Jamison
 673  *
 674  */
 675
 676 /**
 677  * Stage one of clustering.
 678  * @param points float[][] the input points in LAB space
 679  * @param depth int used for recursion, start with 0
 680  * @param clusters ArrayList an Arraylist to store the clusters
 681  * @param limits float[] the cluster diameters
 682  */
 683 static void stageone(std::vector<CLAB> &points,
 684                      int depth,
 685                      std::vector< std::vector<CLAB> > &clusters,
 686                      float *limits)
 687 {
 688     if (points.size() < 1)
 689         return;
 690
 691     int dims=3;
 692     int curdim=depth%dims;
 693     float min = 0.0f;
 694     float max = 0.0f;
 695     if (curdim == 0)
 696         {
 697         min=points[0].C;
 698         max=points[0].C;
 699         // find maximum and minimum
 700         for (unsigned int i=1; i<points.size(); i++) {
 701             if (min>points[i].C)
 702                 min=points[i].C;
 703             if (max<points[i].C)
 704                 max=points[i].C;
 705             }
 706         }
 707     else if (curdim == 1)
 708         {
 709         min=points[0].L;
 710         max=points[0].L;
 711         // find maximum and minimum
 712         for (unsigned int i=1; i<points.size(); i++) {
 713             if (min>points[i].L)
 714                 min=points[i].L;
 715             if (max<points[i].L)
 716                 max=points[i].L;
 717             }
 718         }
 719     else if (curdim == 2)
 720         {
 721         min=points[0].A;
 722         max=points[0].A;
 723         // find maximum and minimum
 724         for (unsigned int i=1; i<points.size(); i++) {
 725             if (min>points[i].A)
 726                 min=points[i].A;
 727             if (max<points[i].A)
 728                 max=points[i].A;
 729             }
 730         }
 731     else if (curdim == 3)
 732         {
 733         min=points[0].B;
 734         max=points[0].B;
 735         // find maximum and minimum
 736         for (unsigned int i=1; i<points.size(); i++) {
 737             if (min>points[i].B)
 738                 min=points[i].B;
 739             if (max<points[i].B)
 740                 max=points[i].B;
 741             }
 742         }
 743
 744     if (max-min>limits[curdim]) { // Split according to Rubner-Rule
 745         // split
 746         float pivotvalue=((max-min)/2.0f)+min;
 747
 748         std::vector<CLAB> smallerpoints; // allocate mem
 749         std::vector<CLAB> biggerpoints;
 750
 751         for (unsigned int i=0; i<points.size(); i++) { // do actual split
 752             float v = 0.0f;
 753             if (curdim==0)
 754                 v = points[i].C;
 755             else if (curdim==1)
 756                 v = points[i].L;
 757             else if (curdim==2)
 758                 v = points[i].A;
 759             else if (curdim==3)
 760                 v = points[i].B;
 761             if (v <= pivotvalue) {
 762                 smallerpoints.push_back(points[i]);
 763             } else {
 764                 biggerpoints.push_back(points[i]);
 765             }
 766         } // create subtrees
 767         stageone(smallerpoints, depth+1, clusters, limits);
 768         stageone(biggerpoints,  depth+1, clusters, limits);
 769     } else { // create leave
 770         clusters.push_back(points);
 771     }
 772 }
 773
 774 /**
 775  * Stage two of clustering.
 776  * @param points float[][] the input points in LAB space
 777  * @param depth int used for recursion, start with 0
 778  * @param clusters ArrayList an Arraylist to store the clusters
 779  * @param limits float[] the cluster diameters
 780  * @param total int the total number of points as given to stageone
 781  * @param threshold should be 0.01 - abstraction threshold
 782  */
 783 static void stagetwo(std::vector<CLAB> &points,
 784                      int depth,
 785                      std::vector< std::vector<CLAB> > &clusters,
 786                      float *limits, int total, float threshold)
 787 {
 788     if (points.size() < 1)
 789         return;
 790
 791     int curdim=depth%3; // without cardinality
 792     float min = 0.0f;
 793     float max = 0.0f;
 794     if (curdim == 0)
 795         {
 796         min=points[0].L;
 797         max=points[0].L;
 798         // find maximum and minimum
 799         for (unsigned int i=1; i<points.size(); i++) {
 800             if (min>points[i].L)
 801                 min=points[i].L;
 802             if (max<points[i].L)
 803                 max=points[i].L;
 804             }
 805         }
 806     else if (curdim == 1)
 807         {
 808         min=points[0].A;
 809         max=points[0].A;
 810         // find maximum and minimum
 811         for (unsigned int i=1; i<points.size(); i++) {
 812             if (min>points[i].A)
 813                 min=points[i].A;
 814             if (max<points[i].A)
 815                 max=points[i].A;
 816             }
 817         }
 818     else if (curdim == 2)
 819         {
 820         min=points[0].B;
 821         max=points[0].B;
 822         // find maximum and minimum
 823         for (unsigned int i=1; i<points.size(); i++) {
 824             if (min>points[i].B)
 825                 min=points[i].B;
 826             if (max<points[i].B)
 827                 max=points[i].B;
 828             }
 829         }
 830
 831
 832     if (max-min>limits[curdim]) { // Split according to Rubner-Rule
 833         // split
 834         float pivotvalue=((max-min)/2.0f)+min;
 835
 836         std::vector<CLAB> smallerpoints; // allocate mem
 837         std::vector<CLAB> biggerpoints;
 838
 839         for (unsigned int i=0; i<points.size(); i++) { // do actual split
 840             float v = 0.0f;
 841             if (curdim==0)
 842                 v = points[i].L;
 843             else if (curdim==1)
 844                 v = points[i].A;
 845             else if (curdim==2)
 846                 v = points[i].B;
 847
 848             if (v <= pivotvalue) {
 849                 smallerpoints.push_back(points[i]);
 850             } else {
 851                 biggerpoints.push_back(points[i]);
 852             }
 853         } // create subtrees
 854         stagetwo(smallerpoints, depth+1, clusters, limits, total, threshold);
 855         stagetwo(biggerpoints,  depth+1, clusters, limits, total, threshold);
 856     } else { // create leave
 857         float sum=0.0;
 858         for (unsigned int i=0; i<points.size(); i++) {
 859             sum+=points[i].B;
 860         }
 861         if (((sum*100.0)/total)>=threshold) {
 862             CLAB point;
 863             for (unsigned int i=0; i<points.size(); i++) {
 864                 point.C += points[i].C;
 865                 point.L += points[i].L;
 866                 point.A += points[i].A;
 867                 point.B += points[i].B;
 868             }
 869             point.C /= points.size();
 870             point.L /= points.size();
 871             point.A /= points.size();
 872             point.B /= points.size();
 873
 874             std::vector<CLAB> newCluster;
 875             newCluster.push_back(point);
 876             clusters.push_back(newCluster);
 877         }
 878     }
 879 }
 880
 881 /**
 882  * Create a color signature for the given set of pixels.
 883  * @param input float[][] a set of pixels in LAB space
 884  * @param length int the number of pixels that should be processed from the input
 885  * @param limits float[] the cluster diameters for each dimension
 886  * @param threshold float the abstraction threshold
 887  * @return float[][] a color siganture containing cluster centroids in LAB space
 888  */
 889 static std::vector<CLAB> createSignature(std::vector<CLAB> &input,
 890                                          float *limits, float threshold)
 891 {
 892     std::vector< std::vector<CLAB> > clusters1;
 893     std::vector< std::vector<CLAB> > clusters2;
 894
 895     stageone(input, 0, clusters1, limits);
 896
 897     std::vector<CLAB> centroids;
 898     for (unsigned int i=0; i<clusters1.size(); i++) {
 899         std::vector<CLAB> cluster = clusters1[i];
 900         CLAB centroid; // +1 for the cardinality
 901         for (unsigned int k=0; k<cluster.size(); k++) {
 902             centroid.L += cluster[k].L;
 903             centroid.A += cluster[k].A;
 904             centroid.B += cluster[k].B;
 905         }
 906         centroid.C =  cluster.size();
 907         centroid.L /= cluster.size();
 908         centroid.A /= cluster.size();
 909         centroid.B /= cluster.size();
 910
 911         centroids.push_back(centroid);
 912     }
 913     stagetwo(centroids, 0, clusters2, limits, input.size(), threshold); // 0.1 -> see paper by tomasi
 914
 915     std::vector<CLAB> res;
 916     for (unsigned int i=0 ; i<clusters2.size() ; i++)
 917         for (unsigned int j=0 ; j<clusters2[i].size() ; j++)
 918             res.push_back(clusters2[i][j]);
 919
 920     return res;
 921 }
 922
 923
 924
 925 //########################################################################
 926 //##  S I O X  S E G M E N T A T O R    (originally SioxSegmentator.java)
 927 //########################################################################
 928
 929 //### NOTE: Doxygen comments are in siox.h
 930
 931 /** Confidence corresponding to a certain foreground region (equals one). */
 932 const float SioxSegmentator::CERTAIN_FOREGROUND_CONFIDENCE=1.0f;
 933
 934 /** Confidence for a region likely being foreground.*/
 935 const float SioxSegmentator::FOREGROUND_CONFIDENCE=0.8f;
 936
 937 /** Confidence for foreground or background type being equally likely.*/
 938 const float SioxSegmentator::UNKNOWN_REGION_CONFIDENCE=0.5f;
 939
 940 /** Confidence for a region likely being background.*/
 941 const float SioxSegmentator::BACKGROUND_CONFIDENCE=0.1f;
 942
 943 /** Confidence corresponding to a certain background reagion (equals zero). */
 944 const float SioxSegmentator::CERTAIN_BACKGROUND_CONFIDENCE=0.0f;
 945
 946
 947 SioxSegmentator::SioxSegmentator(int w, int h, float *limitsArg, int limitsSize)
 948 {
 949
 950     imgWidth   = w;
 951     imgHeight  = h;
 952     labelField = new int[imgWidth*imgHeight];
 953
 954     if (!limitsArg) {
 955         limits = new float[3];
 956         limits[0] = 0.64f;
 957         limits[1] = 1.28f;
 958         limits[2] = 2.56f;
 959     } else {
 960         limits = new float[limitsSize];
 961         for (int i=0 ; i<limitsSize ; i++)
 962             limits[i] = limitsArg[i];
 963     }
 964
 965     float negLimits[3];
 966     negLimits[0] = -limits[0];
 967     negLimits[1] = -limits[1];
 968     negLimits[2] = -limits[2];
 969     clusterSize = sqrEuclidianDist(limits, 3, negLimits);
 970
 971     segmentated=false;
 972
 973     origImage = NULL;
 974 }
 975
 976 SioxSegmentator::~SioxSegmentator()
 977 {
 978     delete labelField;
 979     delete limits;
 980     delete origImage;
 981 }
 982
 983
 984 /**
 985  * Error logging
 986  */
 987 void SioxSegmentator::error(char *fmt, ...)
 988 {
 989     va_list args;
 990     fprintf(stderr, "SioxSegmentator error:");
 991     va_start(args, fmt);
 992     vfprintf(stderr, fmt, args);
 993     va_end(args) ;
 994     fprintf(stderr, "\n");
 995 }
 996
 997 /**
 998  * Trace logging
 999  */
1000 void SioxSegmentator::trace(char *fmt, ...)
1001 {
1002     va_list args;
1003     fprintf(stderr, "SioxSegmentator:");
1004     va_start(args, fmt);
1005     vfprintf(stderr, fmt, args);
1006     va_end(args) ;
1007     fprintf(stderr, "\n");
1008 }
1009
1010
1011
1012 bool SioxSegmentator::segmentate(unsigned long *image, int imageSize,
1013                                  float *cm, int cmSize,
1014                                  int smoothness, double sizeFactorToKeep)
1015 {
1016     segmentated=false;
1017
1018     hs.clear();
1019
1020     // save image for drb
1021     origImage=new long[imageSize];
1022     for (int i=0 ; i<imageSize ; i++)
1023         origImage[i] = image[i];
1024
1025     // create color signatures
1026     for (int i=0; i<cmSize; i++) {
1027         if (cm[i]<=BACKGROUND_CONFIDENCE)
1028             knownBg.push_back(rgbToClab(image[i]));
1029         else if (cm[i]>=FOREGROUND_CONFIDENCE)
1030             knownFg.push_back(rgbToClab(image[i]));
1031     }
1032
1033     bgSignature = createSignature(knownBg, limits, BACKGROUND_CONFIDENCE);
1034     fgSignature = createSignature(knownFg, limits, BACKGROUND_CONFIDENCE);
1035
1036     if (bgSignature.size() < 1) {
1037         // segmentation impossible
1038         return false;
1039     }
1040
1041     // classify using color signatures,
1042     // classification cached in hashmap for drb and speedup purposes
1043     for (int i=0; i<cmSize; i++) {
1044         if (cm[i]>=FOREGROUND_CONFIDENCE) {
1045             cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1046             continue;
1047         }
1048         if (cm[i]>BACKGROUND_CONFIDENCE) {
1049             bool isBackground=true;
1050             std::map<unsigned long, Tupel>::iterator iter = hs.find(i);
1051             Tupel tupel(0.0f, 0, 0.0f, 0);
1052             if (iter == hs.end()) {
1053                 CLAB lab = rgbToClab(image[i]);
1054                 float minBg = sqrEuclidianDist(lab, bgSignature[0]);
1055                 int minIndex=0;
1056                 for (unsigned int j=1; j<bgSignature.size(); j++) {
1057                     float d = sqrEuclidianDist(lab, bgSignature[j]);
1058                     if (d<minBg) {
1059                         minBg=d;
1060                         minIndex=j;
1061                     }
1062                 }
1063                 tupel.minBgDist=minBg;
1064                 tupel.indexMinBg=minIndex;
1065                 float minFg = 1.0e6f;
1066                 minIndex=-1;
1067                 for (unsigned int j=0 ; j<fgSignature.size() ; j++) {
1068                     float d = sqrEuclidianDist(lab, fgSignature[j]);
1069                     if (d<minFg) {
1070                         minFg=d;
1071                         minIndex=j;
1072                     }
1073                 }
1074                 tupel.minFgDist=minFg;
1075                 tupel.indexMinFg=minIndex;
1076                 if (fgSignature.size()==0) {
1077                     isBackground=(minBg<=clusterSize);
1078                     // remove next line to force behaviour of old algorithm
1079                     error("foreground signature does not exist");
1080                     return false;
1081                 } else {
1082                     isBackground=minBg<minFg;
1083                 }
1084                 hs[image[i]] = tupel;
1085             } else {
1086                 isBackground=tupel.minBgDist<=tupel.minFgDist;
1087             }
1088             if (isBackground) {
1089                 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1090             } else {
1091                 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1092             }
1093         } else {
1094             cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1095         }
1096     }
1097
1098     // postprocessing
1099     smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
1100     normalizeMatrix(cm, cmSize);
1101     erode(cm, imgWidth, imgHeight);
1102     keepOnlyLargeComponents(cm, cmSize, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
1103
1104     for (int i=0; i<smoothness; i++) {
1105         smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
1106     }
1107
1108     normalizeMatrix(cm, cmSize);
1109
1110     for (int i=0; i<cmSize; i++) {
1111         if (cm[i]>=UNKNOWN_REGION_CONFIDENCE) {
1112             cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1113         } else {
1114             cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1115         }
1116     }
1117
1118     keepOnlyLargeComponents(cm, cmSize, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
1119     fillColorRegions(cm, cmSize, image);
1120     dilate(cm, imgWidth, imgHeight);
1121
1122     segmentated=true;
1123     return true;
1124 }
1125
1126
1127
1128 void SioxSegmentator::keepOnlyLargeComponents(float *cm, int cmSize,
1129                                               float threshold,
1130                                               double sizeFactorToKeep)
1131 {
1132     int idx = 0;
1133     for (int i=0 ; i<imgHeight ; i++)
1134         for (int j=0 ; j<imgWidth ; j++)
1135             labelField[idx++] = -1;
1136
1137     int curlabel = 0;
1138     int maxregion= 0;
1139     int maxblob  = 0;
1140
1141     // slow but easy to understand:
1142     std::vector<int> labelSizes;
1143     for (int i=0 ; i<cmSize ; i++) {
1144         regionCount=0;
1145         if (labelField[i]==-1 && cm[i]>=threshold) {
1146             regionCount=depthFirstSearch(cm, i, threshold, curlabel++);
1147             labelSizes.push_back(regionCount);
1148         }
1149
1150         if (regionCount>maxregion) {
1151             maxregion=regionCount;
1152             maxblob=curlabel-1;
1153         }
1154     }
1155
1156     for (int i=0 ; i<cmSize ; i++) {
1157         if (labelField[i]!=-1) {
1158             // remove if the component is to small
1159             if (labelSizes[labelField[i]]*sizeFactorToKeep < maxregion)
1160                 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1161
1162             // add maxblob always to foreground
1163             if (labelField[i]==maxblob)
1164                 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1165         }
1166     }
1167 }
1168
1169
1170
1171 int SioxSegmentator::depthFirstSearch(float *cm, int i, float threshold, int curLabel)
1172 {
1173     // stores positions of labeled pixels, where the neighbours
1174     // should still be checked for processing:
1175     std::vector<int> pixelsToVisit;
1176     int componentSize=0;
1177     if (labelField[i]==-1 && cm[i]>=threshold) { // label #i
1178         labelField[i] = curLabel;
1179         ++componentSize;
1180         pixelsToVisit.push_back(i);
1181     }
1182     while (pixelsToVisit.size() > 0) {
1183         int pos=pixelsToVisit[pixelsToVisit.size()-1];
1184         pixelsToVisit.erase(pixelsToVisit.end()-1);
1185         int x=pos%imgWidth;
1186         int y=pos/imgWidth;
1187         // check all four neighbours
1188         int left = pos-1;
1189         if (x-1>=0 && labelField[left]==-1 && cm[left]>=threshold) {
1190             labelField[left]=curLabel;
1191             ++componentSize;
1192             pixelsToVisit.push_back(left);
1193         }
1194         int right = pos+1;
1195         if (x+1<imgWidth && labelField[right]==-1 && cm[right]>=threshold) {
1196             labelField[right]=curLabel;
1197             ++componentSize;
1198             pixelsToVisit.push_back(right);
1199         }
1200         int top = pos-imgWidth;
1201         if (y-1>=0 && labelField[top]==-1 && cm[top]>=threshold) {
1202             labelField[top]=curLabel;
1203             ++componentSize;
1204             pixelsToVisit.push_back(top);
1205         }
1206         int bottom = pos+imgWidth;
1207         if (y+1<imgHeight && labelField[bottom]==-1
1208             && cm[bottom]>=threshold) {
1209             labelField[bottom]=curLabel;
1210             ++componentSize;
1211             pixelsToVisit.push_back(bottom);
1212         }
1213     }
1214     return componentSize;
1215 }
1216
1217 void SioxSegmentator::subpixelRefine(int x, int y, int brushmode,
1218                              float threshold, float *cf, int brushsize)
1219 {
1220     subpixelRefine(x-brushsize, y-brushsize,
1221                    2*brushsize, 2*brushsize,
1222                    brushmode, threshold, cf);
1223 }
1224
1225
1226 bool SioxSegmentator::subpixelRefine(int xa, int ya, int dx, int dy,
1227                                      int brushmode,
1228                                      float threshold, float *cf)
1229 {
1230     if (!segmentated) {
1231         error("no segmentation yet");
1232         return false;
1233     }
1234
1235     int x0 = (xa > 0) ? xa : 0;
1236     int y0 = (ya > 0) ? ya : 0;
1237
1238     int xTo = (imgWidth  - 1 < xa+dx ) ? imgWidth-1  : xa+dx;
1239     int yTo = (imgHeight - 1 < ya+dy ) ? imgHeight-1 : ya+dy;
1240
1241     for (int ey=y0; ey<yTo; ++ey) {
1242         for (int ex=x0; ex<xTo; ++ex) {
1243             /*  we are using a rect, not necessary
1244             if (!area.contains(ex, ey)) {
1245                 continue;
1246             }
1247             */
1248             unsigned long val=origImage[ey*imgWidth+ex];
1249             unsigned long orig=val;
1250             float minDistBg = 0.0f;
1251             float minDistFg = 0.0f;
1252             std::map<unsigned long, Tupel>::iterator iter = hs.find(val);
1253             if (iter != hs.end()) {
1254                 minDistBg=(float) sqrt((float)iter->second.minBgDist);
1255                 minDistFg=(float) sqrt((float)iter->second.minFgDist);
1256             } else {
1257                 continue;
1258             }
1259             if (ADD_EDGE == brushmode) { // handle adder
1260                 if (cf[ey*imgWidth+ex]<FOREGROUND_CONFIDENCE) { // postprocessing wins
1261                   float alpha;
1262                   if (minDistFg==0) {
1263                       alpha=CERTAIN_FOREGROUND_CONFIDENCE;
1264                   } else {
1265                       alpha = (minDistBg/minDistFg < CERTAIN_FOREGROUND_CONFIDENCE) ?
1266                           minDistBg/minDistFg : CERTAIN_FOREGROUND_CONFIDENCE;
1267                   }
1268                   if (alpha<threshold) { // background with certain confidence decided by user.
1269                       alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1270                   }
1271                   val = setAlpha(alpha, orig);
1272                   cf[ey*imgWidth+ex]=alpha;
1273                 }
1274             } else if (SUB_EDGE == brushmode) { // handle substractor
1275                 if (cf[ey*imgWidth+ex]>FOREGROUND_CONFIDENCE) {
1276                     // foreground, we want to take something away
1277                     float alpha;
1278                     if (minDistBg==0) {
1279                         alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1280                     } else {
1281                         alpha=CERTAIN_FOREGROUND_CONFIDENCE-
1282                             (minDistFg/minDistBg < CERTAIN_FOREGROUND_CONFIDENCE) ? // more background -> >1
1283                              minDistFg/minDistBg : CERTAIN_FOREGROUND_CONFIDENCE;
1284                         // bg = gf -> 1
1285                         // more fg -> <1
1286                     }
1287                     if (alpha<threshold) { // background with certain confidence decided by user
1288                         alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1289                     }
1290                     val = setAlpha(alpha, orig);
1291                     cf[ey*imgWidth+ex]=alpha;
1292                 }
1293             } else {
1294                 error("unknown brush mode: %d", brushmode);
1295                 return false;
1296             }
1297         }
1298     }
1299     return true;
1300 }
1301
1302
1303
1304 void SioxSegmentator::fillColorRegions(float *cm, int cmSize, unsigned long *image)
1305 {
1306     int idx = 0;
1307     for (int i=0 ; i<imgHeight ; i++)
1308         for (int j=0 ; i<imgWidth ; j++)
1309             labelField[idx++] = -1;
1310
1311     //int maxRegion=0; // unused now
1312     std::vector<int> pixelsToVisit;
1313     for (int i=0; i<cmSize; i++) { // for all pixels
1314         if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE) {
1315             continue; // already visited or bg
1316         }
1317         int origColor=image[i];
1318         int curLabel=i+1;
1319         labelField[i]=curLabel;
1320         cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1321         // int componentSize = 1;
1322         pixelsToVisit.push_back(i);
1323         // depth first search to fill region
1324         while (pixelsToVisit.size() > 0) {
1325             int pos=pixelsToVisit[pixelsToVisit.size()-1];
1326             pixelsToVisit.erase(pixelsToVisit.end()-1);
1327             int x=pos%imgWidth;
1328             int y=pos/imgWidth;
1329             // check all four neighbours
1330             int left = pos-1;
1331             if (x-1>=0 && labelField[left]==-1
1332                 && labcolordiff(image[left], origColor)<1.0) {
1333                 labelField[left]=curLabel;
1334                 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1335                 // ++componentSize;
1336                 pixelsToVisit.push_back(left);
1337             }
1338             int right = pos+1;
1339             if (x+1<imgWidth && labelField[right]==-1
1340                 && labcolordiff(image[right], origColor)<1.0) {
1341                 labelField[right]=curLabel;
1342                 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1343                 // ++componentSize;
1344                 pixelsToVisit.push_back(right);
1345             }
1346             int top = pos-imgWidth;
1347             if (y-1>=0 && labelField[top]==-1
1348                 && labcolordiff(image[top], origColor)<1.0) {
1349                 labelField[top]=curLabel;
1350                 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1351                 // ++componentSize;
1352                 pixelsToVisit.push_back(top);
1353             }
1354             int bottom = pos+imgWidth;
1355             if (y+1<imgHeight && labelField[bottom]==-1
1356                 && labcolordiff(image[bottom], origColor)<1.0) {
1357                 labelField[bottom]=curLabel;
1358                 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1359                 // ++componentSize;
1360                 pixelsToVisit.push_back(bottom);
1361             }
1362         }
1363         //if (componentSize>maxRegion) {
1364         //    maxRegion=componentSize;
1365         //}
1366     }
1367 }
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382 } //namespace siox
1383 } //namespace org
1384
1385
1386