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     pixelCount = imgWidth * imgHeight;
 953
 954     labelField = new int[pixelCount];
 955
 956     if (!limitsArg) {
 957         limits = new float[3];
 958         limits[0] = 0.64f;
 959         limits[1] = 1.28f;
 960         limits[2] = 2.56f;
 961     } else {
 962         limits = new float[limitsSize];
 963         for (int i=0 ; i<limitsSize ; i++)
 964             limits[i] = limitsArg[i];
 965     }
 966
 967     float negLimits[3];
 968     negLimits[0] = -limits[0];
 969     negLimits[1] = -limits[1];
 970     negLimits[2] = -limits[2];
 971     clusterSize = sqrEuclidianDist(limits, 3, negLimits);
 972
 973     segmentated=false;
 974
 975     origImage = NULL;
 976 }
 977
 978 SioxSegmentator::~SioxSegmentator()
 979 {
 980     delete labelField;
 981     delete limits;
 982     delete origImage;
 983 }
 984
 985
 986 /**
 987  * Error logging
 988  */
 989 void SioxSegmentator::error(char *fmt, ...)
 990 {
 991     va_list args;
 992     fprintf(stderr, "SioxSegmentator error:");
 993     va_start(args, fmt);
 994     vfprintf(stderr, fmt, args);
 995     va_end(args) ;
 996     fprintf(stderr, "\n");
 997 }
 998
 999 /**
1000  * Trace logging
1001  */
1002 void SioxSegmentator::trace(char *fmt, ...)
1003 {
1004     va_list args;
1005     fprintf(stderr, "SioxSegmentator:");
1006     va_start(args, fmt);
1007     vfprintf(stderr, fmt, args);
1008     va_end(args) ;
1009     fprintf(stderr, "\n");
1010 }
1011
1012
1013
1014 bool SioxSegmentator::segmentate(unsigned long *image,
1015                                  float *cm,
1016                                  int smoothness, double sizeFactorToKeep)
1017 {
1018     segmentated=false;
1019
1020     hs.clear();
1021
1022     // save image for drb
1023     origImage=new long[pixelCount];
1024     for (int i=0 ; i<pixelCount ; i++)
1025         origImage[i] = image[i];
1026
1027     // create color signatures
1028     for (int i=0; i<pixelCount; i++) {
1029         if (cm[i]<=BACKGROUND_CONFIDENCE)
1030             knownBg.push_back(rgbToClab(image[i]));
1031         else if (cm[i]>=FOREGROUND_CONFIDENCE)
1032             knownFg.push_back(rgbToClab(image[i]));
1033     }
1034
1035     bgSignature = createSignature(knownBg, limits, BACKGROUND_CONFIDENCE);
1036     fgSignature = createSignature(knownFg, limits, BACKGROUND_CONFIDENCE);
1037
1038     if (bgSignature.size() < 1) {
1039         // segmentation impossible
1040         return false;
1041     }
1042
1043     // classify using color signatures,
1044     // classification cached in hashmap for drb and speedup purposes
1045     for (int i=0; i<pixelCount; i++) {
1046         if (cm[i]>=FOREGROUND_CONFIDENCE) {
1047             cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1048             continue;
1049         }
1050         if (cm[i]>BACKGROUND_CONFIDENCE) {
1051             bool isBackground=true;
1052             std::map<unsigned long, Tupel>::iterator iter = hs.find(i);
1053             Tupel tupel(0.0f, 0, 0.0f, 0);
1054             if (iter == hs.end()) {
1055                 CLAB lab = rgbToClab(image[i]);
1056                 float minBg = sqrEuclidianDist(lab, bgSignature[0]);
1057                 int minIndex=0;
1058                 for (unsigned int j=1; j<bgSignature.size(); j++) {
1059                     float d = sqrEuclidianDist(lab, bgSignature[j]);
1060                     if (d<minBg) {
1061                         minBg=d;
1062                         minIndex=j;
1063                     }
1064                 }
1065                 tupel.minBgDist=minBg;
1066                 tupel.indexMinBg=minIndex;
1067                 float minFg = 1.0e6f;
1068                 minIndex=-1;
1069                 for (unsigned int j=0 ; j<fgSignature.size() ; j++) {
1070                     float d = sqrEuclidianDist(lab, fgSignature[j]);
1071                     if (d<minFg) {
1072                         minFg=d;
1073                         minIndex=j;
1074                     }
1075                 }
1076                 tupel.minFgDist=minFg;
1077                 tupel.indexMinFg=minIndex;
1078                 if (fgSignature.size()==0) {
1079                     isBackground=(minBg<=clusterSize);
1080                     // remove next line to force behaviour of old algorithm
1081                     error("foreground signature does not exist");
1082                     return false;
1083                 } else {
1084                     isBackground=minBg<minFg;
1085                 }
1086                 hs[image[i]] = tupel;
1087             } else {
1088                 isBackground=tupel.minBgDist<=tupel.minFgDist;
1089             }
1090             if (isBackground) {
1091                 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1092             } else {
1093                 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1094             }
1095         } else {
1096             cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1097         }
1098     }
1099
1100     // postprocessing
1101     smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
1102     normalizeMatrix(cm, pixelCount);
1103     erode(cm, imgWidth, imgHeight);
1104     keepOnlyLargeComponents(cm, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
1105
1106     for (int i=0; i<smoothness; i++) {
1107         smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
1108     }
1109
1110     normalizeMatrix(cm, pixelCount);
1111
1112     for (int i=0; i<pixelCount; i++) {
1113         if (cm[i]>=UNKNOWN_REGION_CONFIDENCE) {
1114             cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1115         } else {
1116             cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1117         }
1118     }
1119
1120     keepOnlyLargeComponents(cm, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
1121     fillColorRegions(cm, image);
1122     dilate(cm, imgWidth, imgHeight);
1123
1124     segmentated=true;
1125     return true;
1126 }
1127
1128
1129
1130 void SioxSegmentator::keepOnlyLargeComponents(float *cm,
1131                                               float threshold,
1132                                               double sizeFactorToKeep)
1133 {
1134     int idx = 0;
1135     for (int i=0 ; i<imgHeight ; i++)
1136         for (int j=0 ; j<imgWidth ; j++)
1137             labelField[idx++] = -1;
1138
1139     int curlabel = 0;
1140     int maxregion= 0;
1141     int maxblob  = 0;
1142
1143     // slow but easy to understand:
1144     std::vector<int> labelSizes;
1145     for (int i=0 ; i<pixelCount ; i++) {
1146         regionCount=0;
1147         if (labelField[i]==-1 && cm[i]>=threshold) {
1148             regionCount=depthFirstSearch(cm, i, threshold, curlabel++);
1149             labelSizes.push_back(regionCount);
1150         }
1151
1152         if (regionCount>maxregion) {
1153             maxregion=regionCount;
1154             maxblob=curlabel-1;
1155         }
1156     }
1157
1158     for (int i=0 ; i<pixelCount ; i++) {
1159         if (labelField[i]!=-1) {
1160             // remove if the component is to small
1161             if (labelSizes[labelField[i]]*sizeFactorToKeep < maxregion)
1162                 cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
1163
1164             // add maxblob always to foreground
1165             if (labelField[i]==maxblob)
1166                 cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1167         }
1168     }
1169 }
1170
1171
1172
1173 int SioxSegmentator::depthFirstSearch(float *cm, int i, float threshold, int curLabel)
1174 {
1175     // stores positions of labeled pixels, where the neighbours
1176     // should still be checked for processing:
1177     std::vector<int> pixelsToVisit;
1178     int componentSize=0;
1179     if (labelField[i]==-1 && cm[i]>=threshold) { // label #i
1180         labelField[i] = curLabel;
1181         ++componentSize;
1182         pixelsToVisit.push_back(i);
1183     }
1184     while (pixelsToVisit.size() > 0) {
1185         int pos=pixelsToVisit[pixelsToVisit.size()-1];
1186         pixelsToVisit.erase(pixelsToVisit.end()-1);
1187         int x=pos%imgWidth;
1188         int y=pos/imgWidth;
1189         // check all four neighbours
1190         int left = pos-1;
1191         if (x-1>=0 && labelField[left]==-1 && cm[left]>=threshold) {
1192             labelField[left]=curLabel;
1193             ++componentSize;
1194             pixelsToVisit.push_back(left);
1195         }
1196         int right = pos+1;
1197         if (x+1<imgWidth && labelField[right]==-1 && cm[right]>=threshold) {
1198             labelField[right]=curLabel;
1199             ++componentSize;
1200             pixelsToVisit.push_back(right);
1201         }
1202         int top = pos-imgWidth;
1203         if (y-1>=0 && labelField[top]==-1 && cm[top]>=threshold) {
1204             labelField[top]=curLabel;
1205             ++componentSize;
1206             pixelsToVisit.push_back(top);
1207         }
1208         int bottom = pos+imgWidth;
1209         if (y+1<imgHeight && labelField[bottom]==-1
1210             && cm[bottom]>=threshold) {
1211             labelField[bottom]=curLabel;
1212             ++componentSize;
1213             pixelsToVisit.push_back(bottom);
1214         }
1215     }
1216     return componentSize;
1217 }
1218
1219 void SioxSegmentator::subpixelRefine(int x, int y, int brushmode,
1220                              float threshold, float *cf, int brushsize)
1221 {
1222     subpixelRefine(x-brushsize, y-brushsize,
1223                    2*brushsize, 2*brushsize,
1224                    brushmode, threshold, cf);
1225 }
1226
1227
1228 bool SioxSegmentator::subpixelRefine(int xa, int ya, int dx, int dy,
1229                                      int brushmode,
1230                                      float threshold, float *cf)
1231 {
1232     if (!segmentated) {
1233         error("no segmentation yet");
1234         return false;
1235     }
1236
1237     int x0 = (xa > 0) ? xa : 0;
1238     int y0 = (ya > 0) ? ya : 0;
1239
1240     int xTo = (imgWidth  - 1 < xa+dx ) ? imgWidth-1  : xa+dx;
1241     int yTo = (imgHeight - 1 < ya+dy ) ? imgHeight-1 : ya+dy;
1242
1243     for (int ey=y0; ey<yTo; ++ey) {
1244         for (int ex=x0; ex<xTo; ++ex) {
1245             /*  we are using a rect, not necessary
1246             if (!area.contains(ex, ey)) {
1247                 continue;
1248             }
1249             */
1250             unsigned long val=origImage[ey*imgWidth+ex];
1251             unsigned long orig=val;
1252             float minDistBg = 0.0f;
1253             float minDistFg = 0.0f;
1254             std::map<unsigned long, Tupel>::iterator iter = hs.find(val);
1255             if (iter != hs.end()) {
1256                 minDistBg=(float) sqrt((float)iter->second.minBgDist);
1257                 minDistFg=(float) sqrt((float)iter->second.minFgDist);
1258             } else {
1259                 continue;
1260             }
1261             if (ADD_EDGE == brushmode) { // handle adder
1262                 if (cf[ey*imgWidth+ex]<FOREGROUND_CONFIDENCE) { // postprocessing wins
1263                   float alpha;
1264                   if (minDistFg==0) {
1265                       alpha=CERTAIN_FOREGROUND_CONFIDENCE;
1266                   } else {
1267                       alpha = (minDistBg/minDistFg < CERTAIN_FOREGROUND_CONFIDENCE) ?
1268                           minDistBg/minDistFg : CERTAIN_FOREGROUND_CONFIDENCE;
1269                   }
1270                   if (alpha<threshold) { // background with certain confidence decided by user.
1271                       alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1272                   }
1273                   val = setAlpha(alpha, orig);
1274                   cf[ey*imgWidth+ex]=alpha;
1275                 }
1276             } else if (SUB_EDGE == brushmode) { // handle substractor
1277                 if (cf[ey*imgWidth+ex]>FOREGROUND_CONFIDENCE) {
1278                     // foreground, we want to take something away
1279                     float alpha;
1280                     if (minDistBg==0) {
1281                         alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1282                     } else {
1283                         alpha=CERTAIN_FOREGROUND_CONFIDENCE-
1284                             (minDistFg/minDistBg < CERTAIN_FOREGROUND_CONFIDENCE) ? // more background -> >1
1285                              minDistFg/minDistBg : CERTAIN_FOREGROUND_CONFIDENCE;
1286                         // bg = gf -> 1
1287                         // more fg -> <1
1288                     }
1289                     if (alpha<threshold) { // background with certain confidence decided by user
1290                         alpha=CERTAIN_BACKGROUND_CONFIDENCE;
1291                     }
1292                     val = setAlpha(alpha, orig);
1293                     cf[ey*imgWidth+ex]=alpha;
1294                 }
1295             } else {
1296                 error("unknown brush mode: %d", brushmode);
1297                 return false;
1298             }
1299         }
1300     }
1301     return true;
1302 }
1303
1304
1305
1306 void SioxSegmentator::fillColorRegions(float *cm, unsigned long *image)
1307 {
1308     int idx = 0;
1309     for (int i=0 ; i<imgHeight ; i++)
1310         for (int j=0 ; i<imgWidth ; j++)
1311             labelField[idx++] = -1;
1312
1313     //int maxRegion=0; // unused now
1314     std::vector<int> pixelsToVisit;
1315     for (int i=0; i<pixelCount; i++) { // for all pixels
1316         if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE) {
1317             continue; // already visited or bg
1318         }
1319         int origColor=image[i];
1320         int curLabel=i+1;
1321         labelField[i]=curLabel;
1322         cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
1323         // int componentSize = 1;
1324         pixelsToVisit.push_back(i);
1325         // depth first search to fill region
1326         while (pixelsToVisit.size() > 0) {
1327             int pos=pixelsToVisit[pixelsToVisit.size()-1];
1328             pixelsToVisit.erase(pixelsToVisit.end()-1);
1329             int x=pos%imgWidth;
1330             int y=pos/imgWidth;
1331             // check all four neighbours
1332             int left = pos-1;
1333             if (x-1>=0 && labelField[left]==-1
1334                 && labcolordiff(image[left], origColor)<1.0) {
1335                 labelField[left]=curLabel;
1336                 cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
1337                 // ++componentSize;
1338                 pixelsToVisit.push_back(left);
1339             }
1340             int right = pos+1;
1341             if (x+1<imgWidth && labelField[right]==-1
1342                 && labcolordiff(image[right], origColor)<1.0) {
1343                 labelField[right]=curLabel;
1344                 cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
1345                 // ++componentSize;
1346                 pixelsToVisit.push_back(right);
1347             }
1348             int top = pos-imgWidth;
1349             if (y-1>=0 && labelField[top]==-1
1350                 && labcolordiff(image[top], origColor)<1.0) {
1351                 labelField[top]=curLabel;
1352                 cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
1353                 // ++componentSize;
1354                 pixelsToVisit.push_back(top);
1355             }
1356             int bottom = pos+imgWidth;
1357             if (y+1<imgHeight && labelField[bottom]==-1
1358                 && labcolordiff(image[bottom], origColor)<1.0) {
1359                 labelField[bottom]=curLabel;
1360                 cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
1361                 // ++componentSize;
1362                 pixelsToVisit.push_back(bottom);
1363             }
1364         }
1365         //if (componentSize>maxRegion) {
1366         //    maxRegion=componentSize;
1367         //}
1368     }
1369 }
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384 } //namespace siox
1385 } //namespace org
1386
1387
1388