From: ishmal <ishmal@users.sourceforge.net>
Date: Fri, 24 Mar 2006 14:40:42 +0000 (+0000)
Subject: Add siox
X-Git-Url: https://git.tokkee.org/?a=commitdiff_plain;h=3823b0c996be544aa088e6f03c6e189b66d79a5c;p=inkscape.git

Add siox
---

diff --git a/src/trace/Makefile_insert b/src/trace/Makefile_insert
index 92686be0b..6a2100207 100644
--- a/src/trace/Makefile_insert
+++ b/src/trace/Makefile_insert
@@ -12,8 +12,10 @@ trace_libtrace_a_SOURCES =	            \
 	trace/imagemap-gdk.h	            \
 	trace/imagemap.cpp	            \
 	trace/imagemap.h	            \
-	trace/filterset.cpp	            \
 	trace/filterset.h                   \
+	trace/filterset.cpp	            \
+	trace/siox-segmentator.h	    \
+	trace/siox-segmentator.cpp	    \
         trace/potrace/auxiliary.h           \
         trace/potrace/bitmap.h              \
         trace/potrace/curve.cpp             \
diff --git a/src/trace/siox-segmentator.cpp b/src/trace/siox-segmentator.cpp
new file mode 100644
index 000000000..31d99ca03
--- /dev/null
+++ b/src/trace/siox-segmentator.cpp
@@ -0,0 +1,1371 @@
+/*
+   Copyright 2005, 2006 by Gerald Friedland, Kristian Jantz and Lars Knipping
+
+   Conversion to C++ for Inkscape by Bob Jamison
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+   http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+ */
+#include "siox-segmentator.h"
+
+#include <string>
+
+#include <stdarg.h> //for error() and trace()
+#include <math.h>   //sqrt(), pow(), round(), etc
+
+
+namespace org
+{
+namespace siox
+{
+
+//########################################################################
+//##  U T I L S    (originally Utils.java)
+//########################################################################
+
+/**
+ * Collection of auxiliary image processing methods used by the
+ * SioxSegmentator mainly for postprocessing.
+ *
+ * @author G. Friedland, K. Jantz, L. Knipping
+ * @version 1.05
+ *
+ * Conversion to C++ by Bob Jamison
+ *
+ */
+
+
+/** Caches color conversion values to speed up RGB->CIELAB conversion.*/
+static std::map<long, CLAB> RGB_TO_LAB;
+
+//forward decls
+static void premultiplyMatrix(float alpha, float *cm, int cmSize);
+//static float colordiffsq(long rgb0, long rgb1);
+//static int getAlpha(long argb);
+static int getRed(long rgb);
+static int getGreen(long rgb);
+static int getBlue(long rgb);
+//static long packPixel(int a, int r, int g, int b);
+static CLAB rgbToClab(long rgb);
+
+/**
+ * Applies the morphological dilate operator.
+ *
+ * Can be used to close small holes in the given confidence matrix.
+ *
+ * @param cm Confidence matrix to be processed.
+ * @param xres Horizontal resolution of the matrix.
+ * @param yres Vertical resolution of the matrix.
+ */
+static void dilate(float *cm, int xres, int yres)
+{
+    for (int y=0; y<yres; y++) {
+        for (int x=0; x<xres-1; x++) {
+            int idx=(y*xres)+x;
+            if (cm[idx+1]>cm[idx])
+                cm[idx]=cm[idx+1];
+            }
+        }
+    for (int y=0; y<yres; y++) {
+            for (int x=xres-1; x>=1; x--) {
+            int idx=(y*xres)+x;
+            if (cm[idx-1]>cm[idx])
+                cm[idx]=cm[idx-1];
+        }
+    }
+    for (int y=0; y<yres-1; y++) {
+        for (int x=0; x<xres; x++) {
+        int idx=(y*xres)+x;
+        if (cm[((y+1)*xres)+x] > cm[idx])
+            cm[idx]=cm[((y+1)*xres)+x];
+        }
+    }
+    for (int y=yres-1; y>=1; y--) {
+        for (int x=0; x<xres; x++) {
+        int idx=(y*xres)+x;
+        if (cm[((y-1)*xres)+x] > cm[idx])
+            cm[idx]=cm[((y-1)*xres)+x];
+        }
+    }
+}
+
+/**
+ * Applies the morphological erode operator.
+ *
+ * @param cm Confidence matrix to be processed.
+ * @param xres Horizontal resolution of the matrix.
+ * @param yres Vertical resolution of the matrix.
+ */
+static void erode(float *cm, int xres, int yres)
+{
+    for (int y=0; y<yres; y++) {
+        for (int x=0; x<xres-1; x++) {
+            int idx=(y*xres)+x;
+            if (cm[idx+1] < cm[idx])
+                cm[idx]=cm[idx+1];
+        }
+    }
+    for (int y=0; y<yres; y++) {
+        for (int x=xres-1; x>=1; x--) {
+            int idx=(y*xres)+x;
+            if (cm[idx-1] < cm[idx])
+                cm[idx]=cm[idx-1];
+        }
+    }
+    for (int y=0; y<yres-1; y++) {
+        for (int x=0; x<xres; x++) {
+            int idx=(y*xres)+x;
+            if (cm[((y+1)*xres)+x] < cm[idx])
+                cm[idx]=cm[((y+1)*xres)+x];
+        }
+    }
+    for (int y=yres-1; y>=1; y--) {
+        for (int x=0; x<xres; x++) {
+            int idx=(y*xres)+x;
+            if (cm[((y-1)*xres)+x] < cm[idx])
+                cm[idx]=cm[((y-1)*xres)+x];
+        }
+    }
+}
+
+/**
+ * Normalizes the matrix to values to [0..1].
+ *
+ * @param cm The matrix to be normalized.
+ */
+static void normalizeMatrix(float *cm, int cmSize)
+{
+    float max=0.0f;
+    for (int i=0; i<cmSize; i++) {
+        if (max<cm[i])
+            max=cm[i];
+    }
+    if (max<=0.0)
+        return;
+    else if (max==1.00)
+        return;
+
+    float alpha=1.00f/max;
+    premultiplyMatrix(alpha, cm, cmSize);
+}
+
+/**
+ * Multiplies matrix with the given scalar.
+ *
+ * @param alpha The scalar value.
+ * @param cm The matrix of values be multiplied with alpha.
+ * @param cmSize The matrix length.
+ */
+static void premultiplyMatrix(float alpha, float *cm, int cmSize)
+{
+    for (int i=0; i<cmSize; i++)
+        cm[i]=alpha*cm[i];
+}
+
+/**
+ * Blurs confidence matrix with a given symmetrically weighted kernel.
+ * <P>
+ * In the standard case confidence matrix entries are between 0...1 and
+ * the weight factors sum up to 1.
+ *
+ * @param cm The matrix to be smoothed.
+ * @param xres Horizontal resolution of the matrix.
+ * @param yres Vertical resolution of the matrix.
+ * @param f1 Weight factor for the first pixel.
+ * @param f2 Weight factor for the mid-pixel.
+ * @param f3 Weight factor for the last pixel.
+ */
+static void smoothcm(float *cm, int xres, int yres,
+                     float f1, float f2, float f3)
+{
+    for (int y=0; y<yres; y++) {
+        for (int x=0; x<xres-2; x++) {
+            int idx=(y*xres)+x;
+            cm[idx]=f1*cm[idx]+f2*cm[idx+1]+f3*cm[idx+2];
+        }
+    }
+    for (int y=0; y<yres; y++) {
+        for (int x=xres-1; x>=2; x--) {
+            int idx=(y*xres)+x;
+            cm[idx]=f3*cm[idx-2]+f2*cm[idx-1]+f1*cm[idx];
+        }
+    }
+    for (int y=0; y<yres-2; y++) {
+        for (int x=0; x<xres; x++) {
+            int idx=(y*xres)+x;
+            cm[idx]=f1*cm[idx]+f2*cm[((y+1)*xres)+x]+f3*cm[((y+2)*xres)+x];
+        }
+    }
+    for (int y=yres-1; y>=2; y--) {
+        for (int x=0; x<xres; x++) {
+            int idx=(y*xres)+x;
+            cm[idx]=f3*cm[((y-2)*xres)+x]+f2*cm[((y-1)*xres)+x]+f1*cm[idx];
+        }
+    }
+}
+
+/**
+ * Squared Euclidian distance of p and q.
+ * <P>
+ * Usage hint: When only comparisons between Euclidian distances without
+ * actual values are needed, the squared distance can be preferred
+ * for being faster to compute.
+ *
+ * @param p First euclidian point coordinates.
+ * @param pSize Length of coordinate array.
+ * @param q Second euclidian point coordinates.
+ *        Dimension must not be smaller than that of p.
+ *        Any extra dimensions will be ignored.
+ * @return Squared euclidian distance of p and q.
+ * @see #euclid
+ */
+static float sqrEuclidianDist(float *p, int pSize, float *q)
+{
+    float sum=0;
+    for (int i=0; i<pSize; i++)
+        sum+=(p[i]-q[i])*(p[i]-q[i]);
+    return sum;
+}
+
+/**
+ * Squared Euclidian distance of p and q.
+ * <P>
+ * Usage hint: When only comparisons between Euclidian distances without
+ * actual values are needed, the squared distance can be preferred
+ * for being faster to compute.
+ *
+ * @param p CLAB value
+ * @param q second CLAB value
+ * @return Squared euclidian distance of p and q.
+ * @see #euclid
+ */
+static float sqrEuclidianDist(const CLAB &p, const CLAB &q)
+{
+    float sum=0;
+    sum += (p.C - q.C) * (p.C - q.C);
+    sum += (p.L - q.L) * (p.L - q.L);
+    sum += (p.A - q.A) * (p.A - q.A);
+    sum += (p.B - q.B) * (p.B - q.B);
+    return sum;
+}
+
+/**
+ * Euclidian distance of p and q.
+ *
+ * @param p First euclidian point coordinates.
+ * @param pSize Length of coordinate array.
+ * @param q Second euclidian point coordinates.
+ *        Dimension must not be smaller than that of p.
+ *        Any extra dimensions will be ignored.
+ * @return Squared euclidian distance of p and q.
+ * @see #sqrEuclidianDist
+ */
+/*
+static float euclid(float *p, int pSize, float *q)
+{
+    return (float)sqrt(sqrEuclidianDist(p, pSize, q));
+}
+*/
+
+/**
+ * Computes Euclidian distance of two RGB color values.
+ *
+ * @param rgb0 First color value.
+ * @param rgb1 Second color value.
+ * @return Euclidian distance between the two color values.
+ */
+/*
+static float colordiff(long rgb0, long rgb1)
+{
+    return (float)sqrt(colordiffsq(rgb0, rgb1));
+}
+*/
+
+/**
+ * Computes squared euclidian distance of two RGB color values.
+ * <P>
+ * Note: Faster to compute than colordiff
+ *
+ * @param rgb0 First color value.
+ * @param rgb1 Second color value.
+ * @return Squared Euclidian distance between the two color values.
+ */
+/*
+static float colordiffsq(long rgb0, long rgb1)
+{
+    int rDist=getRed(rgb0)   - getRed(rgb1);
+    int gDist=getGreen(rgb0) - getGreen(rgb1);
+    int bDist=getBlue(rgb0)  - getBlue(rgb1);
+
+    return (float)(rDist*rDist+gDist*gDist+bDist*bDist);
+}
+*/
+
+/**
+ * Averages two ARGB colors.
+ *
+ * @param argb0 First color value.
+ * @param argb1 Second color value.
+ * @return The averaged ARGB color.
+ */
+/*
+static long average(long argb0, long argb1)
+{
+    long ret = packPixel(
+         (getAlpha(argb0) + getAlpha(argb1))/2,
+         (getRed(argb0)   + getRed(argb1)  )/2,
+         (getGreen(argb0) + getGreen(argb1))/2,
+         (getBlue(argb0)  + getBlue(argb1) )/2);
+    return ret;
+}
+*/
+
+/**
+ * Computes squared euclidian distance in CLAB space for two colors
+ * given as RGB values.
+ *
+ * @param rgb0 First color value.
+ * @param rgb1 Second color value.
+ * @return Squared Euclidian distance in CLAB space.
+ */
+static float labcolordiffsq(long rgb1, long rgb2)
+{
+    CLAB c1 = rgbToClab(rgb1);
+    CLAB c2 = rgbToClab(rgb2);
+    float euclid=0.0f;
+    euclid += (c1.L - c2.L) * (c1.L - c2.L);
+    euclid += (c1.A - c2.A) * (c1.A - c2.A);
+    euclid += (c1.B - c2.B) * (c1.B - c2.B);
+    return euclid;
+}
+
+
+/**
+ * Computes squared euclidian distance in CLAB space for two colors
+ * given as RGB values.
+ *
+ * @param rgb0 First color value.
+ * @param rgb1 Second color value.
+ * @return Euclidian distance in CLAB space.
+ */
+static float labcolordiff(long rgb0, long rgb1)
+{
+    return (float)sqrt(labcolordiffsq(rgb0, rgb1));
+}
+
+
+/**
+ * Converts 24-bit RGB values to {l,a,b} float values.
+ * <P>
+ * The conversion used is decribed at
+ * <a href="http://www.easyrgb.com/math.php?MATH=M7#text7">CLAB Conversion</a>
+ * for reference white D65. Note that that the conversion is computational
+ * expensive. Result are cached to speed up further conversion calls.
+ *
+ * @param rgb RGB color value,
+ * @return CLAB color value tripel.
+ */
+static CLAB rgbToClab(long rgb)
+{
+    std::map<long, CLAB>::iterator iter = RGB_TO_LAB.find(rgb);
+    if (iter != RGB_TO_LAB.end())
+        {
+        CLAB res = iter->second;
+        return res;
+        }
+
+    int R=getRed(rgb);
+    int G=getGreen(rgb);
+    int B=getBlue(rgb);
+
+    float var_R=(R/255.0f); //R = From 0 to 255
+    float var_G=(G/255.0f); //G = From 0 to 255
+    float var_B=(B/255.0f); //B = From 0 to 255
+
+    if (var_R>0.04045)
+        var_R=(float) pow((var_R+0.055f)/1.055f, 2.4);
+    else
+        var_R=var_R/12.92f;
+
+    if (var_G>0.04045)
+        var_G=(float) pow((var_G+0.055f)/1.055f, 2.4);
+    else
+        var_G=var_G/12.92f;
+
+    if (var_B>0.04045)
+        var_B=(float) pow((var_B+0.055f)/1.055f, 2.4);
+    else
+        var_B=var_B/12.92f;
+
+    var_R=var_R*100.0f;
+    var_G=var_G*100.0f;
+    var_B=var_B*100.0f;
+
+    // Observer. = 2�, Illuminant = D65
+    float X=var_R*0.4124f + var_G*0.3576f + var_B*0.1805f;
+    float Y=var_R*0.2126f + var_G*0.7152f + var_B*0.0722f;
+    float Z=var_R*0.0193f + var_G*0.1192f + var_B*0.9505f;
+
+    float var_X=X/95.047f;
+    float var_Y=Y/100.0f;
+    float var_Z=Z/108.883f;
+
+    if (var_X>0.008856f)
+        var_X=(float) pow(var_X, 0.3333f);
+    else
+        var_X=(7.787f*var_X)+(16.0f/116.0f);
+
+    if (var_Y>0.008856f)
+        var_Y=(float) pow(var_Y, 0.3333f);
+    else
+        var_Y=(7.787f*var_Y)+(16.0f/116.0f);
+
+    if (var_Z>0.008856f)
+        var_Z=(float) pow(var_Z, 0.3333f);
+    else
+        var_Z=(7.787f*var_Z)+(16.0f/116.0f);
+
+    CLAB lab((116.0f*var_Y)-16.0f , 500.0f*(var_X-var_Y), 200.0f*(var_Y-var_Z));
+
+    RGB_TO_LAB[rgb] = lab;
+
+    return lab;
+}
+
+/**
+ * Converts an CLAB value to a RGB color value.
+ * <P>
+ * This is the reverse operation to rgbToClab.
+ * @param clab CLAB value.
+ * @return RGB value.
+ * @see #rgbToClab
+ */
+/*
+static long clabToRGB(const CLAB &clab)
+{
+    float L=clab.L;
+    float a=clab.A;
+    float b=clab.B;
+
+    float var_Y=(L+16.0f)/116.0f;
+    float var_X=a/500.0f+var_Y;
+    float var_Z=var_Y-b/200.0f;
+
+    float var_yPow3=(float)pow(var_Y, 3.0);
+    float var_xPow3=(float)pow(var_X, 3.0);
+    float var_zPow3=(float)pow(var_Z, 3.0);
+
+    if (var_yPow3>0.008856f)
+        var_Y=var_yPow3;
+    else
+        var_Y=(var_Y-16.0f/116.0f)/7.787f;
+
+    if (var_xPow3>0.008856f)
+        var_X=var_xPow3;
+    else
+        var_X=(var_X-16.0f/116.0f)/7.787f;
+
+    if (var_zPow3>0.008856f)
+        var_Z=var_zPow3;
+    else
+        var_Z=(var_Z-16.0f/116.0f)/7.787f;
+
+    float X=95.047f  * var_X; //ref_X= 95.047  Observer=2�, Illuminant=D65
+    float Y=100.0f   * var_Y; //ref_Y=100.000
+    float Z=108.883f * var_Z; //ref_Z=108.883
+
+    var_X=X/100.0f; //X = From 0 to ref_X
+    var_Y=Y/100.0f; //Y = From 0 to ref_Y
+    var_Z=Z/100.0f; //Z = From 0 to ref_Y
+
+    float var_R=(float)(var_X *  3.2406f + var_Y * -1.5372f + var_Z * -0.4986f);
+    float var_G=(float)(var_X * -0.9689f + var_Y *  1.8758f + var_Z *  0.0415f);
+    float var_B=(float)(var_X *  0.0557f + var_Y * -0.2040f + var_Z *  1.0570f);
+
+    if (var_R>0.0031308f)
+        var_R=(float)(1.055f*pow(var_R, (1.0f/2.4f))-0.055f);
+    else
+        var_R=12.92f*var_R;
+
+    if (var_G>0.0031308f)
+        var_G=(float)(1.055f*pow(var_G, (1.0f/2.4f))-0.055f);
+    else
+        var_G=12.92f*var_G;
+
+    if (var_B>0.0031308f)
+        var_B=(float)(1.055f*pow(var_B, (1.0f/2.4f))-0.055f);
+    else
+        var_B=12.92f*var_B;
+
+    int R = (int)lround(var_R*255.0f);
+    int G = (int)lround(var_G*255.0f);
+    int B = (int)lround(var_B*255.0f);
+
+    return packPixel(0xFF, R, G, B);
+}
+*/
+
+/**
+ * Sets the alpha byte of a pixel.
+ *
+ * Constructs alpha to values from 0 to 255.
+ * @param alpha Alpha value from 0 (transparent) to 255 (opaque).
+ * @param rgb The 24bit rgb color to be combined with the alpga value.
+ * @return An ARBG calor value.
+ */
+static long setAlpha(int alpha, long rgb)
+{
+    if (alpha>255)
+        alpha=0;
+    else if (alpha<0)
+        alpha=0;
+    return (alpha<<24)|(rgb&0xFFFFFF);
+}
+
+/**
+ * Sets the alpha byte of a pixel.
+ *
+ * Constricts alpha to values from 0 to 255.
+ * @param alpha Alpha value from 0.0f (transparent) to 1.0f (opaque).
+ * @param rgb The 24bit rgb color to be combined with the alpga value.
+ * @return An ARBG calor value.
+ */
+static long setAlpha(float alpha, long rgb)
+{
+    return setAlpha((int)(255.0f*alpha), rgb);
+}
+
+/**
+ * Limits the values of a,r,g,b to values from 0 to 255 and puts them
+ * together into an 32 bit integer.
+ *
+ * @param a Alpha part, the first byte.
+ * @param r Red part, the second byte.
+ * @param g Green part, the third byte.
+ * @param b Blue part, the fourth byte.
+ * @return A ARBG value packed to an int.
+ */
+/*
+static long packPixel(int a, int r, int g, int b)
+{
+    if (a<0)
+        a=0;
+    else if (a>255)
+        a=255;
+
+    if (r<0)
+        r=0;
+    else if (r>255)
+        r=255;
+
+    if (g<0)
+        g=0;
+    else if (g>255)
+        g=255;
+
+    if (b<0)
+        b=0;
+    else if (b>255)
+        b=255;
+
+    return (a<<24)|(r<<16)|(g<<8)|b;
+}
+*/
+
+/**
+ * Returns the alpha component of an ARGB value.
+ *
+ * @param argb An ARGB color value.
+ * @return The alpha component, ranging from 0 to 255.
+ */
+/*
+static int getAlpha(long argb)
+{
+    return (argb>>24)&0xFF;
+}
+*/
+
+/**
+ * Returns the red component of an (A)RGB value.
+ *
+ * @param rgb An (A)RGB color value.
+ * @return The red component, ranging from 0 to 255.
+ */
+static int getRed(long rgb)
+{
+    return (rgb>>16)&0xFF;
+}
+
+
+/**
+ * Returns the green component of an (A)RGB value.
+ *
+ * @param rgb An (A)RGB color value.
+ * @return The green component, ranging from 0 to 255.
+ */
+static int getGreen(long rgb)
+{
+    return (rgb>>8)&0xFF;
+}
+
+/**
+ * Returns the blue component of an (A)RGB value.
+ *
+ * @param rgb An (A)RGB color value.
+ * @return The blue component, ranging from 0 to 255.
+ */
+static int getBlue(long rgb)
+{
+    return (rgb)&0xFF;
+}
+
+/**
+ * Returns a string representation of a CLAB value.
+ *
+ * @param clab The CLAB value.
+ * @param clabSize Size of the CLAB value.
+ * @return A string representation of the CLAB value.
+ */
+/*
+static std::string clabToString(const CLAB &clab)
+{
+    std::string buff;
+    char nbuf[60];
+    snprintf(nbuf, 59, "%5.3f, %5.3f, %5.3f", clab.L, clab.A, clab.B);
+    buff = nbuf;
+    return buff;
+}
+*/
+
+//########################################################################
+//##  C O L O R    S I G N A T U R E    (originally ColorSignature.java)
+//########################################################################
+
+/**
+ * Representation of a color signature.
+ * <br><br>
+ * This class implements a clustering algorithm based on a modified kd-tree.
+ * The splitting rule is to simply divide the given interval into two equally
+ * sized subintervals.
+ * In the <code>stageone()</code>, approximate clusters are found by building
+ * up such a tree and stopping when an interval at a node has become smaller
+ * than the allowed cluster diameter, which is given by <code>limits</code>.
+ * At this point, clusters may be split in several nodes.<br>
+ * Therefore, in <code>stagetwo()</code>, nodes that belong to several clusters
+ * are recombined by another k-d tree clustering on the prior cluster
+ * centroids. To guarantee a proper level of abstraction, clusters that contain
+ * less than 0.01% of the pixels of the entire sample are removed. Please
+ * refer to the file NOTICE to get links to further documentation.
+ *
+ * @author Gerald Friedland, Lars Knipping
+ * @version 1.02
+ *
+ * Conversion to C++ by Bob Jamison
+ *
+ */
+
+/**
+ * Stage one of clustering.
+ * @param points float[][] the input points in LAB space
+ * @param depth int used for recursion, start with 0
+ * @param clusters ArrayList an Arraylist to store the clusters
+ * @param limits float[] the cluster diameters
+ */
+static void stageone(std::vector<CLAB> &points,
+                     int depth,
+                     std::vector< std::vector<CLAB> > &clusters,
+                     float *limits)
+{
+    if (points.size() < 1)
+        return;
+
+    int dims=3;
+    int curdim=depth%dims;
+    float min = 0.0f;
+    float max = 0.0f;
+    if (curdim == 0)
+        {
+        min=points[0].C;
+        max=points[0].C;
+        // find maximum and minimum
+        for (unsigned int i=1; i<points.size(); i++) {
+            if (min>points[i].C)
+                min=points[i].C;
+            if (max<points[i].C)
+                max=points[i].C;
+            }
+        }
+    else if (curdim == 1)
+        {
+        min=points[0].L;
+        max=points[0].L;
+        // find maximum and minimum
+        for (unsigned int i=1; i<points.size(); i++) {
+            if (min>points[i].L)
+                min=points[i].L;
+            if (max<points[i].L)
+                max=points[i].L;
+            }
+        }
+    else if (curdim == 2)
+        {
+        min=points[0].A;
+        max=points[0].A;
+        // find maximum and minimum
+        for (unsigned int i=1; i<points.size(); i++) {
+            if (min>points[i].A)
+                min=points[i].A;
+            if (max<points[i].A)
+                max=points[i].A;
+            }
+        }
+    else if (curdim == 3)
+        {
+        min=points[0].B;
+        max=points[0].B;
+        // find maximum and minimum
+        for (unsigned int i=1; i<points.size(); i++) {
+            if (min>points[i].B)
+                min=points[i].B;
+            if (max<points[i].B)
+                max=points[i].B;
+            }
+        }
+
+    if (max-min>limits[curdim]) { // Split according to Rubner-Rule
+        // split
+        float pivotvalue=((max-min)/2.0f)+min;
+
+        std::vector<CLAB> smallerpoints; // allocate mem
+        std::vector<CLAB> biggerpoints;
+
+        for (unsigned int i=0; i<points.size(); i++) { // do actual split
+            float v = 0.0f;
+            if (curdim==0)
+                v = points[i].C;
+            else if (curdim==1)
+                v = points[i].L;
+            else if (curdim==2)
+                v = points[i].A;
+            else if (curdim==3)
+                v = points[i].B;
+            if (v <= pivotvalue) {
+                smallerpoints.push_back(points[i]);
+            } else {
+                biggerpoints.push_back(points[i]);
+            }
+        } // create subtrees
+        stageone(smallerpoints, depth+1, clusters, limits);
+        stageone(biggerpoints,  depth+1, clusters, limits);
+    } else { // create leave
+        clusters.push_back(points);
+    }
+}
+
+/**
+ * Stage two of clustering.
+ * @param points float[][] the input points in LAB space
+ * @param depth int used for recursion, start with 0
+ * @param clusters ArrayList an Arraylist to store the clusters
+ * @param limits float[] the cluster diameters
+ * @param total int the total number of points as given to stageone
+ * @param threshold should be 0.01 - abstraction threshold
+ */
+static void stagetwo(std::vector<CLAB> &points,
+                     int depth,
+                     std::vector< std::vector<CLAB> > &clusters,
+                     float *limits, int total, float threshold)
+{
+    if (points.size() < 1)
+        return;
+
+    int curdim=depth%3; // without cardinality
+    float min = 0.0f;
+    float max = 0.0f;
+    if (curdim == 0)
+        {
+        min=points[0].L;
+        max=points[0].L;
+        // find maximum and minimum
+        for (unsigned int i=1; i<points.size(); i++) {
+            if (min>points[i].L)
+                min=points[i].L;
+            if (max<points[i].L)
+                max=points[i].L;
+            }
+        }
+    else if (curdim == 1)
+        {
+        min=points[0].A;
+        max=points[0].A;
+        // find maximum and minimum
+        for (unsigned int i=1; i<points.size(); i++) {
+            if (min>points[i].A)
+                min=points[i].A;
+            if (max<points[i].A)
+                max=points[i].A;
+            }
+        }
+    else if (curdim == 2)
+        {
+        min=points[0].B;
+        max=points[0].B;
+        // find maximum and minimum
+        for (unsigned int i=1; i<points.size(); i++) {
+            if (min>points[i].B)
+                min=points[i].B;
+            if (max<points[i].B)
+                max=points[i].B;
+            }
+        }
+
+
+    if (max-min>limits[curdim]) { // Split according to Rubner-Rule
+        // split
+        float pivotvalue=((max-min)/2.0f)+min;
+
+        std::vector<CLAB> smallerpoints; // allocate mem
+        std::vector<CLAB> biggerpoints;
+
+        for (unsigned int i=0; i<points.size(); i++) { // do actual split
+            float v = 0.0f;
+	    if (curdim==0)
+	        v = points[i].L;
+	    else if (curdim==1)
+	        v = points[i].A;
+	    else if (curdim==2)
+	        v = points[i].B;
+
+            if (v <= pivotvalue) {
+                smallerpoints.push_back(points[i]);
+            } else {
+                biggerpoints.push_back(points[i]);
+            }
+        } // create subtrees
+        stagetwo(smallerpoints, depth+1, clusters, limits, total, threshold);
+        stagetwo(biggerpoints,  depth+1, clusters, limits, total, threshold);
+    } else { // create leave
+        float sum=0.0;
+        for (unsigned int i=0; i<points.size(); i++) {
+            sum+=points[i].B;
+        }
+        if (((sum*100.0)/total)>=threshold) {
+            CLAB point;
+            for (unsigned int i=0; i<points.size(); i++) {
+                point.C += points[i].C;
+                point.L += points[i].L;
+                point.A += points[i].A;
+                point.B += points[i].B;
+            }
+            point.C /= points.size();
+            point.L /= points.size();
+            point.A /= points.size();
+            point.B /= points.size();
+
+            std::vector<CLAB> newCluster;
+            newCluster.push_back(point);
+            clusters.push_back(newCluster);
+        }
+    }
+}
+
+/**
+ * Create a color signature for the given set of pixels.
+ * @param input float[][] a set of pixels in LAB space
+ * @param length int the number of pixels that should be processed from the input
+ * @param limits float[] the cluster diameters for each dimension
+ * @param threshold float the abstraction threshold
+ * @return float[][] a color siganture containing cluster centroids in LAB space
+ */
+static std::vector<CLAB> createSignature(std::vector<CLAB> &input,
+                                         float *limits, float threshold)
+{
+    std::vector< std::vector<CLAB> > clusters1;
+    std::vector< std::vector<CLAB> > clusters2;
+
+    stageone(input, 0, clusters1, limits);
+
+    std::vector<CLAB> centroids;
+    for (unsigned int i=0; i<clusters1.size(); i++) {
+        std::vector<CLAB> cluster = clusters1[i];
+        CLAB centroid; // +1 for the cardinality
+        for (unsigned int k=0; k<cluster.size(); k++) {
+            centroid.L += cluster[k].L;
+            centroid.A += cluster[k].A;
+            centroid.B += cluster[k].B;
+        }
+        centroid.C =  cluster.size();
+        centroid.L /= cluster.size();
+        centroid.A /= cluster.size();
+        centroid.B /= cluster.size();
+
+        centroids.push_back(centroid);
+    }
+    stagetwo(centroids, 0, clusters2, limits, input.size(), threshold); // 0.1 -> see paper by tomasi
+
+    std::vector<CLAB> res;
+    for (unsigned int i=0 ; i<clusters2.size() ; i++)
+        for (unsigned int j=0 ; j<clusters2[i].size() ; j++)
+            res.push_back(clusters2[i][j]);
+
+    return res;
+}
+
+
+
+//########################################################################
+//##  S I O X  S E G M E N T A T O R    (originally SioxSegmentator.java)
+//########################################################################
+
+//### NOTE: Doxygen comments are in siox-segmentator.h
+
+
+SioxSegmentator::SioxSegmentator(int w, int h, float *limitsArg, int limitsSize)
+{
+
+    imgWidth   = w;
+    imgHeight  = h;
+    labelField = new int[imgWidth*imgHeight];
+
+    if (!limitsArg) {
+        limits = new float[3];
+        limits[0] = 0.64f;
+        limits[1] = 1.28f;
+        limits[2] = 2.56f;
+    } else {
+        limits = new float[limitsSize];
+        for (int i=0 ; i<limitsSize ; i++)
+            limits[i] = limitsArg[i];
+    }
+
+    float negLimits[3];
+    negLimits[0] = -limits[0];
+    negLimits[1] = -limits[1];
+    negLimits[2] = -limits[2];
+    clusterSize = sqrEuclidianDist(limits, 3, negLimits);
+
+    segmentated=false;
+
+    origImage = NULL;
+}
+
+SioxSegmentator::~SioxSegmentator()
+{
+    delete labelField;
+    delete limits;
+    delete origImage;
+}
+
+
+/**
+ * Error logging
+ */
+void SioxSegmentator::error(char *fmt, ...)
+{
+    va_list args;
+    fprintf(stderr, "SioxSegmentator error:");
+    va_start(args, fmt);
+    vfprintf(stderr, fmt, args);
+    va_end(args) ;
+    fprintf(stderr, "\n");
+}
+
+/**
+ * Trace logging
+ */
+void SioxSegmentator::trace(char *fmt, ...)
+{
+    va_list args;
+    fprintf(stderr, "SioxSegmentator:");
+    va_start(args, fmt);
+    vfprintf(stderr, fmt, args);
+    va_end(args) ;
+    fprintf(stderr, "\n");
+}
+
+
+
+bool SioxSegmentator::segmentate(long *image, int imageSize,
+                                 float *cm, int cmSize,
+                                 int smoothness, double sizeFactorToKeep)
+{
+    segmentated=false;
+
+    hs.clear();
+
+    // save image for drb
+    origImage=new long[imageSize];
+    for (int i=0 ; i<imageSize ; i++)
+        origImage[i] = image[i];
+
+    // create color signatures
+    for (int i=0; i<cmSize; i++) {
+        if (cm[i]<=BACKGROUND_CONFIDENCE)
+            knownBg.push_back(rgbToClab(image[i]));
+        else if (cm[i]>=FOREGROUND_CONFIDENCE)
+            knownFg.push_back(rgbToClab(image[i]));
+    }
+
+    bgSignature = createSignature(knownBg, limits, BACKGROUND_CONFIDENCE);
+    fgSignature = createSignature(knownFg, limits, BACKGROUND_CONFIDENCE);
+
+    if (bgSignature.size() < 1) {
+        // segmentation impossible
+        return false;
+    }
+
+    // classify using color signatures,
+    // classification cached in hashmap for drb and speedup purposes
+    for (int i=0; i<cmSize; i++) {
+        if (cm[i]>=FOREGROUND_CONFIDENCE) {
+            cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
+            continue;
+        }
+        if (cm[i]>BACKGROUND_CONFIDENCE) {
+            bool isBackground=true;
+            std::map<long, Tupel>::iterator iter = hs.find(i);
+            Tupel tupel(0.0f, 0, 0.0f, 0);
+            if (iter == hs.end()) {
+                CLAB lab = rgbToClab(image[i]);
+                float minBg = sqrEuclidianDist(lab, bgSignature[0]);
+                int minIndex=0;
+                for (unsigned int j=1; j<bgSignature.size(); j++) {
+                    float d = sqrEuclidianDist(lab, bgSignature[j]);
+                    if (d<minBg) {
+                        minBg=d;
+                        minIndex=j;
+                    }
+                }
+                tupel.minBgDist=minBg;
+                tupel.indexMinBg=minIndex;
+                float minFg = 1.0e6f;
+                minIndex=-1;
+                for (unsigned int j=0 ; j<fgSignature.size() ; j++) {
+                    float d = sqrEuclidianDist(lab, fgSignature[j]);
+                    if (d<minFg) {
+                        minFg=d;
+                        minIndex=j;
+                    }
+                }
+                tupel.minFgDist=minFg;
+                tupel.indexMinFg=minIndex;
+                if (fgSignature.size()==0) {
+                    isBackground=(minBg<=clusterSize);
+                    // remove next line to force behaviour of old algorithm
+                    error("foreground signature does not exist");
+                    return false;
+                } else {
+                    isBackground=minBg<minFg;
+                }
+                hs[image[i]] = tupel;
+            } else {
+                isBackground=tupel.minBgDist<=tupel.minFgDist;
+            }
+            if (isBackground) {
+                cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
+            } else {
+                cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
+            }
+        } else {
+            cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
+        }
+    }
+
+    // postprocessing
+    smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
+    normalizeMatrix(cm, cmSize);
+    erode(cm, imgWidth, imgHeight);
+    keepOnlyLargeComponents(cm, cmSize, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
+
+    for (int i=0; i<smoothness; i++) {
+        smoothcm(cm, imgWidth, imgHeight, 0.33f, 0.33f, 0.33f); // average
+    }
+
+    normalizeMatrix(cm, cmSize);
+
+    for (int i=0; i<cmSize; i++) {
+        if (cm[i]>=UNKNOWN_REGION_CONFIDENCE) {
+            cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
+        } else {
+            cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
+        }
+    }
+
+    keepOnlyLargeComponents(cm, cmSize, UNKNOWN_REGION_CONFIDENCE, sizeFactorToKeep);
+    fillColorRegions(cm, cmSize, image);
+    dilate(cm, imgWidth, imgHeight);
+
+    segmentated=true;
+    return true;
+}
+
+
+
+void SioxSegmentator::keepOnlyLargeComponents(float *cm, int cmSize,
+                                              float threshold,
+                                              double sizeFactorToKeep)
+{
+    int idx = 0;
+    for (int i=0 ; i<imgHeight ; i++)
+        for (int j=0 ; j<imgWidth ; j++)
+            labelField[idx++] = -1;
+
+    int curlabel = 0;
+    int maxregion= 0;
+    int maxblob  = 0;
+
+    // slow but easy to understand:
+    std::vector<int> labelSizes;
+    for (int i=0 ; i<cmSize ; i++) {
+        regionCount=0;
+        if (labelField[i]==-1 && cm[i]>=threshold) {
+            regionCount=depthFirstSearch(cm, i, threshold, curlabel++);
+            labelSizes.push_back(regionCount);
+        }
+
+        if (regionCount>maxregion) {
+            maxregion=regionCount;
+            maxblob=curlabel-1;
+        }
+    }
+
+    for (int i=0 ; i<cmSize ; i++) {
+        if (labelField[i]!=-1) {
+            // remove if the component is to small
+            if (labelSizes[labelField[i]]*sizeFactorToKeep < maxregion)
+                cm[i]=CERTAIN_BACKGROUND_CONFIDENCE;
+
+            // add maxblob always to foreground
+            if (labelField[i]==maxblob)
+                cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
+        }
+    }
+}
+
+
+
+int SioxSegmentator::depthFirstSearch(float *cm, int i, float threshold, int curLabel)
+{
+    // stores positions of labeled pixels, where the neighbours
+    // should still be checked for processing:
+    std::vector<int> pixelsToVisit;
+    int componentSize=0;
+    if (labelField[i]==-1 && cm[i]>=threshold) { // label #i
+        labelField[i] = curLabel;
+        ++componentSize;
+        pixelsToVisit.push_back(i);
+    }
+    while (pixelsToVisit.size() > 0) {
+        int pos=pixelsToVisit[pixelsToVisit.size()-1];
+        pixelsToVisit.erase(pixelsToVisit.end()-1);
+        int x=pos%imgWidth;
+        int y=pos/imgWidth;
+        // check all four neighbours
+        int left = pos-1;
+        if (x-1>=0 && labelField[left]==-1 && cm[left]>=threshold) {
+            labelField[left]=curLabel;
+            ++componentSize;
+            pixelsToVisit.push_back(left);
+        }
+        int right = pos+1;
+        if (x+1<imgWidth && labelField[right]==-1 && cm[right]>=threshold) {
+            labelField[right]=curLabel;
+            ++componentSize;
+            pixelsToVisit.push_back(right);
+        }
+        int top = pos-imgWidth;
+        if (y-1>=0 && labelField[top]==-1 && cm[top]>=threshold) {
+            labelField[top]=curLabel;
+            ++componentSize;
+            pixelsToVisit.push_back(top);
+        }
+        int bottom = pos+imgWidth;
+        if (y+1<imgHeight && labelField[bottom]==-1
+            && cm[bottom]>=threshold) {
+            labelField[bottom]=curLabel;
+            ++componentSize;
+            pixelsToVisit.push_back(bottom);
+        }
+    }
+    return componentSize;
+}
+
+void SioxSegmentator::subpixelRefine(int x, int y, int brushmode,
+                             float threshold, float *cf, int brushsize)
+{
+    subpixelRefine(x-brushsize, y-brushsize,
+                   2*brushsize, 2*brushsize,
+                   brushmode, threshold, cf);
+}
+
+
+bool SioxSegmentator::subpixelRefine(int xa, int ya, int dx, int dy,
+                                     int brushmode,
+                                     float threshold, float *cf)
+{
+    if (!segmentated) {
+        error("no segmentation yet");
+        return false;
+    }
+
+    int x0 = (xa > 0) ? xa : 0;
+    int y0 = (ya > 0) ? ya : 0;
+
+    int xTo = (imgWidth  - 1 < xa+dx ) ? imgWidth-1  : xa+dx;
+    int yTo = (imgHeight - 1 < ya+dy ) ? imgHeight-1 : ya+dy;
+
+    for (int ey=y0; ey<yTo; ++ey) {
+        for (int ex=x0; ex<xTo; ++ex) {
+            /*  we are using a rect, not necessary
+            if (!area.contains(ex, ey)) {
+                continue;
+            }
+            */
+            long val=origImage[ey*imgWidth+ex];
+            long orig=val;
+            float minDistBg = 0.0f;
+            float minDistFg = 0.0f;
+            std::map<long, Tupel>::iterator iter = hs.find(val);
+            if (iter != hs.end()) {
+                minDistBg=(float) sqrt((float)iter->second.minBgDist);
+                minDistFg=(float) sqrt((float)iter->second.minFgDist);
+            } else {
+                continue;
+            }
+            if (ADD_EDGE == brushmode) { // handle adder
+                if (cf[ey*imgWidth+ex]<FOREGROUND_CONFIDENCE) { // postprocessing wins
+                  float alpha;
+                  if (minDistFg==0) {
+                      alpha=CERTAIN_FOREGROUND_CONFIDENCE;
+                  } else {
+                      alpha = (minDistBg/minDistFg < CERTAIN_FOREGROUND_CONFIDENCE) ?
+                          minDistBg/minDistFg : CERTAIN_FOREGROUND_CONFIDENCE;
+                  }
+                  if (alpha<threshold) { // background with certain confidence decided by user.
+                      alpha=CERTAIN_BACKGROUND_CONFIDENCE;
+                  }
+                  val = setAlpha(alpha, orig);
+                  cf[ey*imgWidth+ex]=alpha;
+                }
+            } else if (SUB_EDGE == brushmode) { // handle substractor
+                if (cf[ey*imgWidth+ex]>FOREGROUND_CONFIDENCE) {
+                    // foreground, we want to take something away
+                    float alpha;
+                    if (minDistBg==0) {
+                        alpha=CERTAIN_BACKGROUND_CONFIDENCE;
+                    } else {
+                        alpha=CERTAIN_FOREGROUND_CONFIDENCE-
+                            (minDistFg/minDistBg < CERTAIN_FOREGROUND_CONFIDENCE) ? // more background -> >1
+                             minDistFg/minDistBg : CERTAIN_FOREGROUND_CONFIDENCE;
+                        // bg = gf -> 1
+                        // more fg -> <1
+                    }
+                    if (alpha<threshold) { // background with certain confidence decided by user
+                        alpha=CERTAIN_BACKGROUND_CONFIDENCE;
+                    }
+                    val = setAlpha(alpha, orig);
+                    cf[ey*imgWidth+ex]=alpha;
+                }
+            } else {
+                error("unknown brush mode: %d", brushmode);
+                return false;
+            }
+        }
+    }
+    return true;
+}
+
+
+
+void SioxSegmentator::fillColorRegions(float *cm, int cmSize, long *image)
+{
+    int idx = 0;
+    for (int i=0 ; i<imgHeight ; i++)
+        for (int j=0 ; i<imgWidth ; j++)
+            labelField[idx++] = -1;
+
+    //int maxRegion=0; // unused now
+    std::vector<int> pixelsToVisit;
+    for (int i=0; i<cmSize; i++) { // for all pixels
+        if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE) {
+            continue; // already visited or bg
+        }
+        int origColor=image[i];
+        int curLabel=i+1;
+        labelField[i]=curLabel;
+        cm[i]=CERTAIN_FOREGROUND_CONFIDENCE;
+        // int componentSize = 1;
+        pixelsToVisit.push_back(i);
+        // depth first search to fill region
+        while (pixelsToVisit.size() > 0) {
+            int pos=pixelsToVisit[pixelsToVisit.size()-1];
+            pixelsToVisit.erase(pixelsToVisit.end()-1);
+            int x=pos%imgWidth;
+            int y=pos/imgWidth;
+            // check all four neighbours
+            int left = pos-1;
+            if (x-1>=0 && labelField[left]==-1
+                && labcolordiff(image[left], origColor)<1.0) {
+                labelField[left]=curLabel;
+                cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
+                // ++componentSize;
+                pixelsToVisit.push_back(left);
+            }
+            int right = pos+1;
+            if (x+1<imgWidth && labelField[right]==-1
+                && labcolordiff(image[right], origColor)<1.0) {
+                labelField[right]=curLabel;
+                cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
+                // ++componentSize;
+                pixelsToVisit.push_back(right);
+            }
+            int top = pos-imgWidth;
+            if (y-1>=0 && labelField[top]==-1
+                && labcolordiff(image[top], origColor)<1.0) {
+                labelField[top]=curLabel;
+                cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
+                // ++componentSize;
+                pixelsToVisit.push_back(top);
+            }
+            int bottom = pos+imgWidth;
+            if (y+1<imgHeight && labelField[bottom]==-1
+                && labcolordiff(image[bottom], origColor)<1.0) {
+                labelField[bottom]=curLabel;
+                cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
+                // ++componentSize;
+                pixelsToVisit.push_back(bottom);
+            }
+        }
+        //if (componentSize>maxRegion) {
+        //    maxRegion=componentSize;
+        //}
+    }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+} //namespace siox
+} //namespace org
+
+
+
diff --git a/src/trace/siox-segmentator.h b/src/trace/siox-segmentator.h
new file mode 100644
index 000000000..4d92a9182
--- /dev/null
+++ b/src/trace/siox-segmentator.h
@@ -0,0 +1,396 @@
+#ifndef __SIOX_SEGMENTATOR_H__
+#define __SIOX_SEGMENTATOR_H__
+/*
+   Copyright 2005, 2006 by Gerald Friedland, Kristian Jantz and Lars Knipping
+
+   Conversion to C++ for Inkscape by Bob Jamison
+
+   Licensed under the Apache License, Version 2.0 (the "License");
+   you may not use this file except in compliance with the License.
+   You may obtain a copy of the License at
+
+   http://www.apache.org/licenses/LICENSE-2.0
+
+   Unless required by applicable law or agreed to in writing, software
+   distributed under the License is distributed on an "AS IS" BASIS,
+   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+   See the License for the specific language governing permissions and
+   limitations under the License.
+ */
+
+#include <map>
+#include <vector>
+
+namespace org
+{
+namespace siox
+{
+
+/**
+ * Image segmentator based on
+ *<em>SIOX: Simple Interactive Object Extraction</em>.
+ * <P>
+ * To segmentate an image one has to perform the following steps.
+ * <OL><LI>Construct an instance of <code>SioxSegmentator</code>.
+ * </LI><LI>Create a confidence matrix, where each entry marks its
+ *      corresponding image pixel to belong to the foreground, to the
+ *      background, or being of unknown type.
+ * </LI><LI>Call <code>segmentate</code> on the image with the confidence
+ *       matrix. This stores the result as new foreground confidence into
+ *       the confidence matrix, with each entry being either
+ *       zero (<code>CERTAIN_BACKGROUND_CONFIDENCE</code>) or one
+ *       (<code>CERTAIN_FOREGROUND_CONFIDENCE</code>).
+ * </LI><LI>Optionally call <code>subpixelRefine</code> to areas
+ *      where pixels contain both foreground and background (e.g.
+ *      object borders or highly detailed features like flowing hairs).
+ *      The pixel are then assigned confidence values bwetween zero and
+ *      one to give them a measure of "foregroundness".
+ *      This step may be repeated as often as needed.
+ * </LI></OL>
+ * <P>
+ * For algorithm documentation refer to
+ * G. Friedland, K. Jantz, L. Knipping, R. Rojas:<i>
+ * Image Segmentation by Uniform Color Clustering
+ *  -- Approach and Benchmark Results</i>,
+ * <A HREF="http://www.inf.fu-berlin.de/inst/pubs/tr-b-05-07.pdf">Technical Report B-05-07</A>,
+ * Department of Computer Science, Freie Universitaet Berlin, June 2005.<br>
+ * <P>
+ * See <A HREF="http://www.siox.org/" target="_new">http://www.siox.org</A> for more information.<br>
+ * <P>
+ * Algorithm idea by Gerald Friedland.
+ *
+ * @author Gerald Friedland, Kristian Jantz, Lars Knipping
+ * @version 1.12
+ */
+
+/**
+ * Helper class for storing the minimum distances to a cluster centroid
+ * in background and foreground and the index to the centroids in each
+ * signature for a given color.
+ */
+class Tupel {
+public:
+
+    Tupel()
+        {
+        minBgDist  = 0.0f;
+        indexMinBg = 0;
+        minFgDist  = 0.0f;
+        indexMinFg = 0;
+        }
+    Tupel(float minBgDistArg, long indexMinBgArg,
+          float minFgDistArg, long indexMinFgArg)
+        {
+	minBgDist  = minBgDistArg;
+	indexMinBg = indexMinBgArg;
+	minFgDist  = minFgDistArg;
+	indexMinFg = indexMinFgArg;
+        }
+    Tupel(const Tupel &other)
+        {
+	minBgDist  = other.minBgDist;
+	indexMinBg = other.indexMinBg;
+	minFgDist  = other.minFgDist;
+	indexMinFg = other.indexMinFg;
+        }
+    Tupel &operator=(const Tupel &other)
+        {
+	minBgDist  = other.minBgDist;
+	indexMinBg = other.indexMinBg;
+	minFgDist  = other.minFgDist;
+	indexMinFg = other.indexMinFg;
+	return *this;
+        }
+    virtual ~Tupel()
+        {}
+
+    float minBgDist;
+    long  indexMinBg;
+    float minFgDist;
+    long  indexMinFg;
+
+ };
+
+
+class CLAB
+{
+public:
+    CLAB()
+        {
+        C = L = A = B = 0.0f;
+        }
+    CLAB(float lArg, float aArg, float bArg)
+        {
+        C = 0.0f;
+        L = lArg;
+        A = aArg;
+        B = bArg;
+        }
+    CLAB(const CLAB &other)
+        {
+        C = other.C;
+        L = other.L;
+        A = other.A;
+        B = other.B;
+        }
+    CLAB &operator=(const CLAB &other)
+        {
+        C = other.C;
+        L = other.L;
+        A = other.A;
+        B = other.B;
+        return *this;
+        }
+    virtual ~CLAB()
+        {}
+
+    float C;
+    float L;
+    float A;
+    float B;
+};
+
+
+class SioxSegmentator
+{
+public:
+
+    /** Confidence corresponding to a certain foreground region (equals one). */
+    static const float CERTAIN_FOREGROUND_CONFIDENCE=1.0f;
+
+    /** Confidence for a region likely being foreground.*/
+    static const float FOREGROUND_CONFIDENCE=0.8f;
+
+    /** Confidence for foreground or background type being equally likely.*/
+    static const float UNKNOWN_REGION_CONFIDENCE=0.5f;
+
+    /** Confidence for a region likely being background.*/
+    static const float BACKGROUND_CONFIDENCE=0.1f;
+
+    /** Confidence corresponding to a certain background reagion (equals zero). */
+    static const float CERTAIN_BACKGROUND_CONFIDENCE=0.0f;
+
+
+    /**
+     * Constructs a SioxSegmentator Object to be used for image segmentation.
+     *
+     * @param w X resolution of the image to be segmentated.
+     * @param h Y resolution of the image to be segmentated.
+     * @param limits Size of the cluster on LAB axises.
+     *        If <code>null</code>, the default value {0.64f,1.28f,2.56f}
+     *        is used.
+     */
+    SioxSegmentator(int w, int h, float *limitsArg, int limitsSize);
+
+    /**
+     * Destructor
+     */
+    virtual ~SioxSegmentator();
+
+
+    /**
+     * Segmentates the given image with information from the confidence
+     * matrix.
+     * <P>
+     * The confidence entries  of <code>BACKGROUND_CONFIDENCE</code> or less
+     * are mark known background pixel for the segmentation, those
+     * of at least <code>FOREGROUND_CONFIDENCE</code> mark known
+     * foreground pixel for the segmentation. Any other entry is treated
+     * as region of unknown affiliation.
+     * <P>
+     * As result, each pixel is classified either as foregroound or
+     * background, stored back into its <code>cm</code> entry as confidence
+     * <code>CERTAIN_FOREGROUND_CONFIDENCE</code> or
+     * <code>CERTAIN_BACKGROUND_CONFIDENCE</code>.
+     *
+     * @param image Pixel data of the image to be segmentated.
+     *        Every integer represents one ARGB-value.
+     * @param imageSize number of values in image
+     * @param cm Confidence matrix specifying the probability of an image
+     *        belonging to the foreground before and after the segmentation.
+     * @param smoothness Number of smoothing steps in the post processing.
+     * @param sizeFactorToKeep Segmentation retains the largest connected
+     *        foreground component plus any component with size at least
+     *        <CODE>sizeOfLargestComponent/sizeFactorToKeep</CODE>.
+     * @return <CODE>true</CODE> if the segmentation algorithm succeeded,
+     *         <CODE>false</CODE> if segmentation is impossible
+     */
+    bool segmentate(long *image, int imageSize,
+                    float *cm, int cmSize,
+                    int smoothness, double sizeFactorToKeep);
+
+    /**
+     * Clears given confidence matrix except entries for the largest connected
+     * component and every component with
+     * <CODE>size*sizeFactorToKeep >= sizeOfLargestComponent</CODE>.
+     *
+     * @param cm  Confidence matrix to be analysed
+     * @param threshold Pixel visibility threshold.
+     *        Exactly those cm entries larger than threshold are considered
+     *        to be a "visible" foreground pixel.
+     * @param sizeFactorToKeep This method keeps the largest connected
+     *        component plus any component with size at least
+     *        <CODE>sizeOfLargestComponent/sizeFactorToKeep</CODE>.
+     */
+    void keepOnlyLargeComponents(float *cm, int cmSize,
+                                 float threshold,
+                                 double sizeFactorToKeep);
+
+    /**
+     * Depth first search pixels in a foreground component.
+     *
+     * @param cm confidence matrix to be searched.
+     * @param i starting position as index to confidence matrix.
+     * @param threshold defines the minimum value at which a pixel is
+     *        considered foreground.
+     * @param curlabel label no of component.
+     * @return size in pixel of the component found.
+     */
+    int depthFirstSearch(float *cm, int i, float threshold, int curLabel);
+
+    /**
+     * Refines the classification stored in the confidence matrix by modifying
+     * the confidences for regions which have characteristics to both
+     * foreground and background if they fall into the specified square.
+     * <P>
+     * The can be used in displaying the image by assigning the alpha values
+     * of the pixels according to the confidence entries.
+     * <P>
+     * In the algorithm descriptions and examples GUIs this step is referrered
+     * to as <EM>Detail Refinement (Brush)</EM>.
+     *
+     * @param x Horizontal coordinate of the squares center.
+     * @param y Vertical coordinate of the squares center.
+     * @param brushmode Mode of the refinement applied, <CODE>ADD_EDGE</CODE>
+     *        or <CODE>SUB_EDGE</CODE>. Add mode only modifies pixels
+     *        formerly classified as background, sub mode only those
+     *        formerly classified as foreground.
+     * @param threshold Threshold for the add and sub refinement, deciding
+     *        at the confidence level to stop at.
+     * @param cf The confidence matrix to modify, generated by
+     *        <CODE>segmentate</CODE>, possibly already refined by privious
+     *        calls to <CODE>subpixelRefine</CODE>.
+     * @param brushsize Halfed diameter of the square shaped brush.
+     *
+     * @see #segmentate
+     */
+    void SioxSegmentator::subpixelRefine(int x, int y, int brushmode,
+                             float threshold, float *cf, int brushsize);
+
+    /**
+     * Refines the classification stored in the confidence matrix by modifying
+     * the confidences for regions which have characteristics to both
+     * foreground and background if they fall into the specified area.
+     * <P>
+     * The can be used in displaying the image by assigning the alpha values
+     * of the pixels according to the confidence entries.
+     * <P>
+     * In the algorithm descriptions and examples GUIs this step is referrered
+     * to as <EM>Detail Refinement (Brush)</EM>.
+     *
+     * @param area Area in which the reworking of the segmentation is
+     *        applied to.
+     * @param brushmode Mode of the refinement applied, <CODE>ADD_EDGE</CODE>
+     *        or <CODE>SUB_EDGE</CODE>. Add mode only modifies pixels
+     *        formerly classified as background, sub mode only those
+     *        formerly classified as foreground.
+     * @param threshold Threshold for the add and sub refinement, deciding
+     *        at the confidence level to stop at.
+     * @param cf The confidence matrix to modify, generated by
+     *        <CODE>segmentate</CODE>, possibly already refined by privious
+     *        calls to <CODE>subpixelRefine</CODE>.
+     *
+     * @see #segmentate
+     */
+    bool SioxSegmentator::subpixelRefine(int xa, int ya, int dx, int dy,
+                                     int brushmode,
+                                     float threshold, float *cf);
+    /**
+     * A region growing algorithms used to fill up the confidence matrix
+     * with <CODE>CERTAIN_FOREGROUND_CONFIDENCE</CODE> for corresponding
+     * areas of equal colors.
+     * <P>
+     * Basically, the method works like the <EM>Magic Wand<EM> with a
+     * tolerance threshold of zero.
+     *
+     * @param cm confidence matrix to be searched
+     * @param image image to be searched
+     */
+    void fillColorRegions(float *cm, int cmSize, long *image);
+
+private:
+
+    /**
+     * Prevent this from being used
+     */
+    SioxSegmentator();
+
+    /** error logging **/
+    void error(char *format, ...);
+
+    /** trace logging **/
+    void trace(char *format, ...);
+
+    typedef enum
+        {
+        ADD_EDGE, /** Add mode for the subpixel refinement. */
+        SUB_EDGE  /** Subtract mode for the subpixel refinement. */
+        } BrushMode;
+
+    // instance fields:
+
+    /** Horizontal resolution of the image to be segmentated. */
+    int imgWidth;
+
+    /** Vertical resolution of the image to be segmentated. */
+    int imgHeight;
+
+    /** Stores component label (index) by pixel it belongs to. */
+    int *labelField;
+
+    /**
+     * LAB color values of pixels that are definitly known background.
+     * Entries are of form {l,a,b}.
+     */
+    std::vector<CLAB> knownBg;
+
+    /**
+     * LAB color values of pixels that are definitly known foreground.
+     * Entries are of form {l,a,b}.
+     */
+    std::vector<CLAB> knownFg;
+
+    /** Holds background signature (a characteristic subset of the bg.) */
+    std::vector<CLAB> bgSignature;
+
+    /** Holds foreground signature (a characteristic subset of the fg).*/
+    std::vector<CLAB> fgSignature;
+
+    /** Size of cluster on lab axis. */
+    float *limits;
+
+    /** Maximum distance of two lab values. */
+    float clusterSize;
+
+    /**
+     * Stores Tupels for fast access to nearest background/foreground pixels.
+     */
+    std::map<long, Tupel> hs;
+
+    /** Size of the biggest blob.*/
+    int regionCount;
+
+    /** Copy of the original image, needed for detail refinement. */
+    long *origImage;
+    long origImageSize;
+
+    /** A flag that stores if the segmentation algorithm has already ran.*/
+    bool segmentated;
+
+};
+
+} //namespace siox
+} //namespace org
+
+#endif /* __SIOX_SEGMENTATOR_H__ */
+