From cca66b9ec4e494c1d919bff0f71a820d8afab1fa Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sun, 7 Apr 2024 20:24:48 +0200 Subject: Adding upstream version 1.2.2. Signed-off-by: Daniel Baumann --- src/trace/siox.cpp | 1736 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1736 insertions(+) create mode 100644 src/trace/siox.cpp (limited to 'src/trace/siox.cpp') diff --git a/src/trace/siox.cpp b/src/trace/siox.cpp new file mode 100644 index 0000000..49ecfae --- /dev/null +++ b/src/trace/siox.cpp @@ -0,0 +1,1736 @@ +// SPDX-License-Identifier: GPL-2.0-or-later +/* + Copyright 2005, 2006 by Gerald Friedland, Kristian Jantz and Lars Knipping + + Conversion to C++ for Inkscape by Bob Jamison + + Released under GNU GPL v2+, read the file 'COPYING' for more information. + */ +#include "siox.h" + +#include +#include +#include +#include +#include + + +namespace org +{ + +namespace siox +{ + + + +//######################################################################## +//# C L A B +//######################################################################## + +/** + * Convert integer A, R, G, B values into an pixel value. + */ +static unsigned long getRGB(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; +} + + + +/** + * Convert float A, R, G, B values (0.0-1.0) into an pixel value. + */ +static unsigned long getRGB(float a, float r, float g, float b) +{ + return getRGB((int)(a * 256.0), + (int)(r * 256.0), + (int)(g * 256.0), + (int)(b * 256.0)); +} + + + +//######################################### +//# Root approximations for large speedup. +//# By njh! +//######################################### +static const int ROOT_TAB_SIZE = 16; +static float cbrt_table[ROOT_TAB_SIZE +1]; + +double CieLab::cbrt(double x) +{ + double y = cbrt_table[int(x*ROOT_TAB_SIZE )]; // assuming x \in [0, 1] + y = (2.0 * y + x/(y*y))/3.0; + y = (2.0 * y + x/(y*y))/3.0; // polish twice + return y; +} + +static float qn_table[ROOT_TAB_SIZE +1]; + +double CieLab::qnrt(double x) +{ + double y = qn_table[int(x*ROOT_TAB_SIZE )]; // assuming x \in [0, 1] + double Y = y*y; + y = (4.0*y + x/(Y*Y))/5.0; + Y = y*y; + y = (4.0*y + x/(Y*Y))/5.0; // polish twice + return y; +} + +double CieLab::pow24(double x) +{ + double onetwo = x*qnrt(x); + return onetwo*onetwo; +} + + +static bool _clab_inited_ = false; +void CieLab::init() +{ + if (!_clab_inited_) + { + cbrt_table[0] = pow(float(1)/float(ROOT_TAB_SIZE*2), 0.3333); + qn_table[0] = pow(float(1)/float(ROOT_TAB_SIZE*2), 0.2); + for(int i = 1; i < ROOT_TAB_SIZE +1; i++) + { + cbrt_table[i] = pow(float(i)/float(ROOT_TAB_SIZE), 0.3333); + qn_table[i] = pow(float(i)/float(ROOT_TAB_SIZE), 0.2); + } + _clab_inited_ = true; + } +} + + + +/** + * Construct this CieLab from an ARGB value + */ +CieLab::CieLab(unsigned long rgb) +{ + init(); + + int ir = (rgb>>16) & 0xff; + int ig = (rgb>> 8) & 0xff; + int ib = (rgb ) & 0xff; + + float fr = ((float)ir) / 255.0; + float fg = ((float)ig) / 255.0; + float fb = ((float)ib) / 255.0; + + //printf("fr:%f fg:%f fb:%f\n", fr, fg, fb); + if (fr > 0.04045) + //fr = (float) pow((fr + 0.055) / 1.055, 2.4); + fr = (float) pow24((fr + 0.055) / 1.055); + else + fr = fr / 12.92; + + if (fg > 0.04045) + //fg = (float) pow((fg + 0.055) / 1.055, 2.4); + fg = (float) pow24((fg + 0.055) / 1.055); + else + fg = fg / 12.92; + + if (fb > 0.04045) + //fb = (float) pow((fb + 0.055) / 1.055, 2.4); + fb = (float) pow24((fb + 0.055) / 1.055); + else + fb = fb / 12.92; + + // Use white = D65 + const float x = fr * 0.4124 + fg * 0.3576 + fb * 0.1805; + const float y = fr * 0.2126 + fg * 0.7152 + fb * 0.0722; + const float z = fr * 0.0193 + fg * 0.1192 + fb * 0.9505; + + float vx = x / 0.95047; + float vy = y; + float vz = z / 1.08883; + + //printf("vx:%f vy:%f vz:%f\n", vx, vy, vz); + if (vx > 0.008856) + //vx = (float) pow(vx, 0.3333); + vx = (float) cbrt(vx); + else + vx = (7.787 * vx) + (16.0 / 116.0); + + if (vy > 0.008856) + //vy = (float) pow(vy, 0.3333); + vy = (float) cbrt(vy); + else + vy = (7.787 * vy) + (16.0 / 116.0); + + if (vz > 0.008856) + //vz = (float) pow(vz, 0.3333); + vz = (float) cbrt(vz); + else + vz = (7.787 * vz) + (16.0 / 116.0); + + C = 0; + L = 116.0 * vy - 16.0; + A = 500.0 * (vx - vy); + B = 200.0 * (vy - vz); +} + + + +/** + * Return this CieLab's value converted to an ARGB value + */ +unsigned long CieLab::toRGB() +{ + float vy = (L + 16.0) / 116.0; + float vx = A / 500.0 + vy; + float vz = vy - B / 200.0; + + float vx3 = vx * vx * vx; + float vy3 = vy * vy * vy; + float vz3 = vz * vz * vz; + + if (vy3 > 0.008856) + vy = vy3; + else + vy = (vy - 16.0 / 116.0) / 7.787; + + if (vx3 > 0.008856) + vx = vx3; + else + vx = (vx - 16.0 / 116.0) / 7.787; + + if (vz3 > 0.008856) + vz = vz3; + else + vz = (vz - 16.0 / 116.0) / 7.787; + + vx *= 0.95047; //use white = D65 + vz *= 1.08883; + + float vr =(float)(vx * 3.2406 + vy * -1.5372 + vz * -0.4986); + float vg =(float)(vx * -0.9689 + vy * 1.8758 + vz * 0.0415); + float vb =(float)(vx * 0.0557 + vy * -0.2040 + vz * 1.0570); + + if (vr > 0.0031308) + vr = (float)(1.055 * pow(vr, (1.0 / 2.4)) - 0.055); + else + vr = 12.92 * vr; + + if (vg > 0.0031308) + vg = (float)(1.055 * pow(vg, (1.0 / 2.4)) - 0.055); + else + vg = 12.92 * vg; + + if (vb > 0.0031308) + vb = (float)(1.055 * pow(vb, (1.0 / 2.4)) - 0.055); + else + vb = 12.92 * vb; + + return getRGB(0.0, vr, vg, vb); +} + + +/** + * Squared Euclidean distance between this and another color + */ +float CieLab::diffSq(const CieLab &other) +{ + float sum=0.0; + sum += (L - other.L) * (L - other.L); + sum += (A - other.A) * (A - other.A); + sum += (B - other.B) * (B - other.B); + return sum; +} + +/** + * Computes squared euclidean distance in CieLab space for two colors + * given as RGB values. + */ +float CieLab::diffSq(unsigned int rgb1, unsigned int rgb2) +{ + CieLab c1(rgb1); + CieLab c2(rgb2); + float euclid = c1.diffSq(c2); + return euclid; +} + + +/** + * Computes squared euclidean distance in CieLab space for two colors + * given as RGB values. + */ +float CieLab::diff(unsigned int rgb0, unsigned int rgb1) +{ + return (float) sqrt(diffSq(rgb0, rgb1)); +} + + + +//######################################################################## +//# T U P E L +//######################################################################## + +/** + * 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) + = default; + virtual ~Tupel() + = default; + + float minBgDist; + long indexMinBg; + float minFgDist; + long indexMinFg; + + }; + + + +//######################################################################## +//# S I O X I M A G E +//######################################################################## + +/** + * Create an image with the given width and height + */ +SioxImage::SioxImage(unsigned int widthArg, unsigned int heightArg) +{ + init(widthArg, heightArg); +} + +/** + * Copy constructor + */ +SioxImage::SioxImage(const SioxImage &other) +{ + pixdata = nullptr; + cmdata = nullptr; + assign(other); +} + +/** + * Assignment + */ +SioxImage &SioxImage::operator=(const SioxImage &other) +{ + assign(other); + return *this; +} + + +/** + * Clean up after use. + */ +SioxImage::~SioxImage() +{ + if (pixdata) delete[] pixdata; + if (cmdata) delete[] cmdata; +} + +/** + * Error logging + */ +void SioxImage::error(const char *fmt, ...) +{ + char msgbuf[256]; + va_list args; + va_start(args, fmt); + vsnprintf(msgbuf, 255, fmt, args); + va_end(args) ; +#ifdef HAVE_GLIB + g_warning("SioxImage error: %s\n", msgbuf); +#else + fprintf(stderr, "SioxImage error: %s\n", msgbuf); +#endif +} + + +/** + * Returns true if the previous operation on this image + * was successful, else false. + */ +bool SioxImage::isValid() +{ + return valid; +} + +/** + * Sets whether an operation was successful, and whether + * this image should be considered a valid one. + * was successful, else false. + */ +void SioxImage::setValid(bool val) +{ + valid = val; +} + + +/** + * Set a pixel at the x,y coordinates to the given value. + * If the coordinates are out of range, do nothing. + */ +void SioxImage::setPixel(unsigned int x, + unsigned int y, + unsigned int pixval) +{ + if (x >= width || y >= height) + { + error("setPixel: out of bounds (%d,%d)/(%d,%d)", + x, y, width, height); + return; + } + unsigned long offset = width * y + x; + pixdata[offset] = pixval; +} + +/** + * Set a pixel at the x,y coordinates to the given r, g, b values. + * If the coordinates are out of range, do nothing. + */ +void SioxImage::setPixel(unsigned int x, unsigned int y, + unsigned int a, + unsigned int r, + unsigned int g, + unsigned int b) +{ + if (x >= width || y >= height) + { + error("setPixel: out of bounds (%d,%d)/(%d,%d)", + x, y, width, height); + return; + } + unsigned long offset = width * y + x; + unsigned int pixval = ((a << 24) & 0xff000000) | + ((r << 16) & 0x00ff0000) | + ((g << 8) & 0x0000ff00) | + ((b ) & 0x000000ff); + pixdata[offset] = pixval; +} + + + +/** + * Get a pixel at the x,y coordinates given. If + * the coordinates are out of range, return 0; + */ +unsigned int SioxImage::getPixel(unsigned int x, unsigned int y) +{ + if (x >= width || y >= height) + { + error("getPixel: out of bounds (%d,%d)/(%d,%d)", + x, y, width, height); + return 0L; + } + unsigned long offset = width * y + x; + return pixdata[offset]; +} + +/** + * Return the image data buffer + */ +unsigned int *SioxImage::getImageData() +{ + return pixdata; +} + +/** + * Set a confidence value at the x,y coordinates to the given value. + * If the coordinates are out of range, do nothing. + */ +void SioxImage::setConfidence(unsigned int x, + unsigned int y, + float confval) +{ + if (x >= width || y >= height) + { + error("setConfidence: out of bounds (%d,%d)/(%d,%d)", + x, y, width, height); + return; + } + unsigned long offset = width * y + x; + cmdata[offset] = confval; +} + +/** + * Get a confidence valueat the x,y coordinates given. If + * the coordinates are out of range, return 0; + */ +float SioxImage::getConfidence(unsigned int x, unsigned int y) +{ + if (x >= width || y >= height) + { + error("getConfidence: out of bounds (%d,%d)/(%d,%d)", + x, y, width, height); + return 0.0; + } + unsigned long offset = width * y + x; + return cmdata[offset]; +} + +/** + * Return the confidence data buffer + */ +float *SioxImage::getConfidenceData() +{ + return cmdata; +} + + +/** + * Return the width of this image + */ +int SioxImage::getWidth() +{ + return width; +} + +/** + * Return the height of this image + */ +int SioxImage::getHeight() +{ + return height; +} + +/** + * Initialize values. Used by constructors + */ +void SioxImage::init(unsigned int widthArg, unsigned int heightArg) +{ + valid = true; + width = widthArg; + height = heightArg; + imageSize = width * height; + pixdata = new unsigned int[imageSize]; + cmdata = new float[imageSize]; + for (unsigned long i=0 ; i>24) & 0xff; + unsigned int r = ((rgb>>16) & 0xff); + unsigned int g = ((rgb>> 8) & 0xff); + unsigned int b = ((rgb ) & 0xff); + fputc((unsigned char) r, f); + fputc((unsigned char) g, f); + fputc((unsigned char) b, f); + } + } + fclose(f); + return true; +} + + +#ifdef HAVE_GLIB + +/** + * Create an image from a GdkPixbuf + */ +SioxImage::SioxImage(GdkPixbuf *buf) +{ + if (!buf) + return; + + unsigned int width = gdk_pixbuf_get_width(buf); + unsigned int height = gdk_pixbuf_get_height(buf); + init(width, height); //DO THIS NOW!! + + + guchar *pixldata = gdk_pixbuf_get_pixels(buf); + int rowstride = gdk_pixbuf_get_rowstride(buf); + int n_channels = gdk_pixbuf_get_n_channels(buf); + + //### Fill in the cells with RGB values + int row = 0; + for (unsigned int y=0 ; y> 16) & 0xff;//r + p[1] = (rgb >> 8) & 0xff;//g + p[2] = (rgb ) & 0xff;//b + if ( n_channels > 3 ) + { + p[3] = (rgb >> 24) & 0xff;//a + } + p += n_channels; + } + row += rowstride; + } + return buf; +} + +#endif /* GLIB */ + + + + +//######################################################################## +//# S I O X +//######################################################################## + +//############## +//## PUBLIC +//############## + +/** + * Confidence corresponding to a certain foreground region (equals one). + */ +const float Siox::CERTAIN_FOREGROUND_CONFIDENCE=1.0f; + +/** + * Confidence for a region likely being foreground. + */ +const float Siox::FOREGROUND_CONFIDENCE=0.8f; + +/** + * Confidence for foreground or background type being equally likely. + */ +const float Siox::UNKNOWN_REGION_CONFIDENCE=0.5f; + +/** + * Confidence for a region likely being background. + */ +const float Siox::BACKGROUND_CONFIDENCE=0.1f; + +/** + * Confidence corresponding to a certain background reagion (equals zero). + */ +const float Siox::CERTAIN_BACKGROUND_CONFIDENCE=0.0f; + +/** + * Construct a Siox engine + */ +Siox::Siox() : + sioxObserver(nullptr), + keepGoing(true), + width(0), + height(0), + pixelCount(0), + image(nullptr), + cm(nullptr), + labelField(nullptr) +{ + init(); +} + +/** + * Construct a Siox engine + */ +Siox::Siox(SioxObserver *observer) : + sioxObserver(observer), + keepGoing(true), + width(0), + height(0), + pixelCount(0), + image(nullptr), + cm(nullptr), + labelField(nullptr) +{ + init(); +} + + +/** + * + */ +Siox::~Siox() +{ + cleanup(); +} + + +/** + * Error logging + */ +void Siox::error(const char *fmt, ...) +{ + char msgbuf[256]; + va_list args; + va_start(args, fmt); + vsnprintf(msgbuf, 255, fmt, args); + va_end(args) ; +#ifdef HAVE_GLIB + g_warning("Siox error: %s\n", msgbuf); +#else + fprintf(stderr, "Siox error: %s\n", msgbuf); +#endif +} + +/** + * Trace logging + */ +void Siox::trace(const char *fmt, ...) +{ + char msgbuf[256]; + va_list args; + va_start(args, fmt); + vsnprintf(msgbuf, 255, fmt, args); + va_end(args) ; +#ifdef HAVE_GLIB + g_message("Siox: %s\n", msgbuf); +#else + fprintf(stdout, "Siox: %s\n", msgbuf); +#endif +} + + + +/** + * Progress reporting + */ +bool Siox::progressReport(float percentCompleted) +{ + if (!sioxObserver) + return true; + + bool ret = sioxObserver->progress(percentCompleted); + + if (!ret) + { + trace("User selected abort"); + keepGoing = false; + } + + return ret; +} + + + + +/** + * Extract the foreground of the original image, according + * to the values in the confidence matrix. + */ +SioxImage Siox::extractForeground(const SioxImage &originalImage, + unsigned int backgroundFillColor) +{ + trace("### Start"); + + init(); + keepGoing = true; + + SioxImage workImage = originalImage; + + //# fetch some info from the image + width = workImage.getWidth(); + height = workImage.getHeight(); + pixelCount = width * height; + image = workImage.getImageData(); + cm = workImage.getConfidenceData(); + labelField = new int[pixelCount]; + + trace("### Creating signatures"); + + //#### create color signatures + std::vector knownBg; + std::vector knownFg; + CieLab *imageClab = new CieLab[pixelCount]; + for (unsigned long i=0 ; i= FOREGROUND_CONFIDENCE) + knownFg.push_back(lab); + } + + /* + std::vector imageClab; + for (int y = 0 ; y < workImage.getHeight() ; y++) + for (int x = 0 ; x < workImage.getWidth() ; x++) + { + float cm = workImage.getConfidence(x, y); + unsigned int pix = workImage.getPixel(x, y); + CieLab lab(pix); + imageClab.push_back(lab); + if (cm <= BACKGROUND_CONFIDENCE) + knownBg.push_back(lab); //note: uses CieLab(rgb) + else if (cm >= FOREGROUND_CONFIDENCE) + knownFg.push_back(lab); + } + */ + + if (!progressReport(10.0)) + { + error("User aborted"); + workImage.setValid(false); + delete[] imageClab; + delete[] labelField; + return workImage; + } + + trace("knownBg:%u knownFg:%u", static_cast(knownBg.size()), static_cast(knownFg.size())); + + + std::vector bgSignature ; + if (!colorSignature(knownBg, bgSignature, 3)) + { + error("Could not create background signature"); + workImage.setValid(false); + delete[] imageClab; + delete[] labelField; + return workImage; + } + + if (!progressReport(30.0)) + { + error("User aborted"); + workImage.setValid(false); + delete[] imageClab; + delete[] labelField; + return workImage; + } + + + std::vector fgSignature ; + if (!colorSignature(knownFg, fgSignature, 3)) + { + error("Could not create foreground signature"); + workImage.setValid(false); + delete[] imageClab; + delete[] labelField; + return workImage; + } + + //trace("### bgSignature:%d", bgSignature.size()); + + if (bgSignature.empty()) + { + // segmentation impossible + error("Signature size is < 1. Segmentation is impossible"); + workImage.setValid(false); + delete[] imageClab; + delete[] labelField; + return workImage; + } + + if (!progressReport(30.0)) + { + error("User aborted"); + workImage.setValid(false); + delete[] imageClab; + delete[] labelField; + return workImage; + } + + + // classify using color signatures, + // classification cached in hashmap for drb and speedup purposes + trace("### Analyzing image"); + + std::map hs; + + unsigned int progressResolution = pixelCount / 10; + + for (unsigned int i=0; i= FOREGROUND_CONFIDENCE) + { + cm[i] = CERTAIN_FOREGROUND_CONFIDENCE; + } + else if (cm[i] <= BACKGROUND_CONFIDENCE) + { + cm[i] = CERTAIN_BACKGROUND_CONFIDENCE; + } + else // somewhere in between + { + bool isBackground = true; + std::map::iterator iter = hs.find(i); + if (iter != hs.end()) //found + { + Tupel tupel = iter->second; + isBackground = tupel.minBgDist <= tupel.minFgDist; + } + else + { + CieLab lab = imageClab[i]; + float minBg = lab.diffSq(bgSignature[0]); + int minIndex = 0; + for (unsigned int j=1; j= UNKNOWN_REGION_CONFIDENCE) + cm[i] = CERTAIN_FOREGROUND_CONFIDENCE; + else + cm[i] = CERTAIN_BACKGROUND_CONFIDENCE; + } + + keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.5/*sizeFactorToKeep*/); + fillColorRegions(); + dilate(cm, width, height); + + if (!progressReport(100.0)) + { + error("User aborted"); + delete[] labelField; + workImage.setValid(false); + return workImage; + } + + + //#### We are done. Now clear everything but the background + for (unsigned long i = 0; i max) max = curval; + } + + //Do the Rubner-rule split (sounds like a dance) + if (max - min > limits[currentDim]) + { + float pivotPoint = (min + max) / 2.0; //average + unsigned int left = leftBase; + unsigned int right = rightBase - 1; + + //# partition points according to the dimension + while (true) + { + while ( true ) + { + point = points[left]; + if (point(currentDim) > pivotPoint) + break; + left++; + } + while ( true ) + { + point = points[right]; + if (point(currentDim) <= pivotPoint) + break; + right--; + } + + if (left > right) + break; + + point = points[left]; + points[left] = points[right]; + points[right] = point; + + left++; + right--; + } + + //# Recurse and create sub-trees + colorSignatureStage1(points, leftBase, left, + recursionDepth + 1, clusterCount, dims); + colorSignatureStage1(points, left, rightBase, + recursionDepth + 1, clusterCount, dims); + } + else + { + //create a leaf + CieLab newpoint; + + newpoint.C = rightBase - leftBase; + + for (; leftBase < rightBase ; leftBase++) + { + newpoint.add(points[leftBase]); + } + + //printf("clusters:%d\n", *clusters); + + if (newpoint.C != 0) + newpoint.mul(1.0 / (float)newpoint.C); + points[*clusterCount] = newpoint; + (*clusterCount)++; + } +} + + + +/** + * Stage 2 of the color signature work + */ +void Siox::colorSignatureStage2(CieLab *points, + unsigned int leftBase, + unsigned int rightBase, + unsigned int recursionDepth, + unsigned int *clusterCount, + const float threshold, + const unsigned int dims) +{ + unsigned int currentDim = recursionDepth % dims; + CieLab point = points[leftBase]; + float min = point(currentDim); + float max = min; + + for (unsigned int i = leftBase+ 1; i < rightBase; i++) + { + point = points[i]; + float curval = point(currentDim); + if (curval < min) min = curval; + if (curval > max) max = curval; + } + + //Do the Rubner-rule split (sounds like a dance) + if (max - min > limits[currentDim]) + { + float pivotPoint = (min + max) / 2.0; //average + unsigned int left = leftBase; + unsigned int right = rightBase - 1; + + //# partition points according to the dimension + while (true) + { + while ( true ) + { + point = points[left]; + if (point(currentDim) > pivotPoint) + break; + left++; + } + while ( true ) + { + point = points[right]; + if (point(currentDim) <= pivotPoint) + break; + right--; + } + + if (left > right) + break; + + point = points[left]; + points[left] = points[right]; + points[right] = point; + + left++; + right--; + } + + //# Recurse and create sub-trees + colorSignatureStage2(points, leftBase, left, + recursionDepth + 1, clusterCount, threshold, dims); + colorSignatureStage2(points, left, rightBase, + recursionDepth + 1, clusterCount, threshold, dims); + } + else + { + //### Create a leaf + unsigned int sum = 0; + for (unsigned int i = leftBase; i < rightBase; i++) + sum += points[i].C; + + if ((float)sum >= threshold) + { + float scale = (float)(rightBase - leftBase); + CieLab newpoint; + + for (; leftBase < rightBase; leftBase++) + newpoint.add(points[leftBase]); + + if (scale != 0.0) + newpoint.mul(1.0 / scale); + points[*clusterCount] = newpoint; + (*clusterCount)++; + } + } +} + + + +/** + * Main color signature method + */ +bool Siox::colorSignature(const std::vector &inputVec, + std::vector &result, + const unsigned int dims) +{ + + unsigned int length = inputVec.size(); + + if (length < 1) // no error. just don't do anything + return true; + + CieLab *input = new CieLab [length]; + + if (!input) + { + error("Could not allocate buffer for signature"); + return false; + } + for (unsigned int i=0 ; i < length ; i++) + { + input[i] = inputVec[i]; + } + + unsigned int stage1length = 0; + colorSignatureStage1(input, 0, length, 0, &stage1length, dims); + + unsigned int stage2length = 0; + colorSignatureStage2(input, 0, stage1length, 0, &stage2length, length * 0.001, dims); + + result.clear(); + for (unsigned int i=0 ; i < stage2length ; i++) + { + result.push_back(input[i]); + } + + delete[] input; + + return true; +} + + + +/** + * + */ +void Siox::keepOnlyLargeComponents(float threshold, + double sizeFactorToKeep) +{ + for (unsigned long idx = 0 ; idx labelSizes; + for (unsigned long i=0 ; i= threshold) + { + regionCount = depthFirstSearch(i, threshold, curlabel++); + labelSizes.push_back(regionCount); + } + + if (regionCount>maxregion) + { + maxregion = regionCount; + maxblob = curlabel-1; + } + } + + for (unsigned int i=0 ; i pixelsToVisit; + int componentSize = 0; + + if (labelField[startPos]==-1 && cm[startPos]>=threshold) + { + labelField[startPos] = curLabel; + componentSize++; + pixelsToVisit.push_back(startPos); + } + + + while (!pixelsToVisit.empty()) + { + int pos = pixelsToVisit[pixelsToVisit.size() - 1]; + pixelsToVisit.erase(pixelsToVisit.end() - 1); + unsigned int x = pos % width; + unsigned int y = pos / width; + + // check all four neighbours + int left = pos - 1; + if (((int)x)-1>=0 && labelField[left]==-1 && cm[left]>=threshold) + { + labelField[left]=curLabel; + componentSize++; + pixelsToVisit.push_back(left); + } + + int right = pos + 1; + if (x+1 < width && labelField[right]==-1 && cm[right]>=threshold) + { + labelField[right]=curLabel; + componentSize++; + pixelsToVisit.push_back(right); + } + + int top = pos - width; + if (((int)y)-1>=0 && labelField[top]==-1 && cm[top]>=threshold) + { + labelField[top]=curLabel; + componentSize++; + pixelsToVisit.push_back(top); + } + + int bottom = pos + width; + if (y+1 < height && labelField[bottom]==-1 + && cm[bottom]>=threshold) + { + labelField[bottom]=curLabel; + componentSize++; + pixelsToVisit.push_back(bottom); + } + + } + return componentSize; +} + + + +/** + * + */ +void Siox::fillColorRegions() +{ + for (unsigned long idx = 0 ; idx pixelsToVisit; + for (unsigned long i=0; i= 0 && labelField[left] == -1 + && CieLab::diff(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 < width && labelField[right]==-1 + && CieLab::diff(image[right], origColor)<1.0) + { + labelField[right]=curLabel; + cm[right]=CERTAIN_FOREGROUND_CONFIDENCE; + // ++componentSize; + pixelsToVisit.push_back(right); + } + int top = pos - width; + if (((int)y)-1>=0 && labelField[top]==-1 + && CieLab::diff(image[top], origColor)<1.0) + { + labelField[top]=curLabel; + cm[top]=CERTAIN_FOREGROUND_CONFIDENCE; + // ++componentSize; + pixelsToVisit.push_back(top); + } + int bottom = pos + width; + if (y+1 < height && labelField[bottom]==-1 + && CieLab::diff(image[bottom], origColor)<1.0) + { + labelField[bottom]=curLabel; + cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE; + // ++componentSize; + pixelsToVisit.push_back(bottom); + } + } + //if (componentSize>maxRegion) { + // maxRegion=componentSize; + //} + } +} + + + + +/** + * Applies the morphological dilate operator. + * + * Can be used to close small holes in the given confidence matrix. + */ +void Siox::dilate(float *cm, int xres, int yres) +{ + + for (int y=0; ycm[idx]) + cm[idx]=cm[idx+1]; + } + } + + for (int y=0; y=1; x--) + { + int idx=(y*xres)+x; + if (cm[idx-1]>cm[idx]) + cm[idx]=cm[idx-1]; + } + } + + for (int y=0; y cm[idx]) + cm[idx]=cm[((y+1)*xres)+x]; + } + } + + for (int y=yres-1; y>=1; y--) + { + for (int x=0; x cm[idx]) + cm[idx]=cm[((y-1)*xres)+x]; + } + } +} + +/** + * Applies the morphological erode operator. + */ +void Siox::erode(float *cm, int xres, int yres) +{ + for (int y=0; y=1; x--) + { + int idx=(y*xres)+x; + if (cm[idx-1] < cm[idx]) + cm[idx]=cm[idx-1]; + } + } + for (int y=0; y=1; y--) + { + for (int x=0; x max) max=cm[i]; + + //good to use STL, but max() is not iterative + //float max = *std::max(cm, cm + cmSize); + + if (max<=0.0 || max==1.0) + return; + + float alpha=1.00f/max; + premultiplyMatrix(alpha, cm, cmSize); +} + +/** + * Multiplies matrix with the given scalar. + */ +void Siox::premultiplyMatrix(float alpha, float *cm, int cmSize) +{ + for (int i=0; i=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=2; y--) + { + for (int x=0; x