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