diff options
Diffstat (limited to 'src/trace/siox.cpp')
| -rw-r--r-- | src/trace/siox.cpp | 1736 |
1 files changed, 1736 insertions, 0 deletions
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 <cmath> +#include <cstdarg> +#include <map> +#include <algorithm> +#include <cstdlib> + + +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<imageSize ; i++) + { + pixdata[i] = 0; + cmdata[i] = 0.0; + } +} + +/** + * Assign values to that of another + */ +void SioxImage::assign(const SioxImage &other) +{ + if (pixdata) delete[] pixdata; + if (cmdata) delete[] cmdata; + valid = other.valid; + width = other.width; + height = other.height; + imageSize = width * height; + pixdata = new unsigned int[imageSize]; + cmdata = new float[imageSize]; + for (unsigned long i=0 ; i<imageSize ; i++) + { + pixdata[i] = other.pixdata[i]; + cmdata[i] = other.cmdata[i]; + } +} + + +/** + * Write the image to a PPM file + */ +bool SioxImage::writePPM(const std::string &fileName) +{ + + FILE *f = fopen(fileName.c_str(), "wb"); + if (!f) + return false; + + fprintf(f, "P6 %u %u 255\n", width, height); + + for (unsigned int y=0 ; y<height; y++) + { + for (unsigned int x=0 ; x<width ; x++) + { + unsigned int rgb = getPixel(x, y); + //unsigned int alpha = (rgb>>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<height ; y++) + { + guchar *p = pixldata + row; + for (unsigned int x=0 ; x<width ; x++) + { + int r = (int)p[0]; + int g = (int)p[1]; + int b = (int)p[2]; + int alpha = (int)p[3]; + + setPixel(x, y, alpha, r, g, b); + p += n_channels; + } + row += rowstride; + } + +} + + +/** + * Create a GdkPixbuf from this image + */ +GdkPixbuf *SioxImage::getGdkPixbuf() +{ + bool has_alpha = true; + int n_channels = has_alpha ? 4 : 3; + + guchar *pixdata = (guchar *) + malloc(sizeof(guchar) * width * height * n_channels); + if (!pixdata) + { + error("SioxImage::getGdkPixbuf: can not allocate memory for %d x %d x %d image.", width, height, n_channels); + return nullptr; + } + + int rowstride = width * n_channels; + + GdkPixbuf *buf = gdk_pixbuf_new_from_data(pixdata, + GDK_COLORSPACE_RGB, + has_alpha, 8, width, height, + rowstride, (GdkPixbufDestroyNotify)free, nullptr); + + //### Fill in the cells with RGB values + int row = 0; + for (unsigned int y=0 ; y < height ; y++) + { + guchar *p = pixdata + row; + for (unsigned x=0 ; x < width ; x++) + { + unsigned int rgb = getPixel(x, y); + p[0] = (rgb >> 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<CieLab> knownBg; + std::vector<CieLab> knownFg; + CieLab *imageClab = new CieLab[pixelCount]; + for (unsigned long i=0 ; i<pixelCount ; i++) + { + float conf = cm[i]; + unsigned int pix = image[i]; + CieLab lab(pix); + imageClab[i] = lab; + if (conf <= BACKGROUND_CONFIDENCE) + knownBg.push_back(lab); + else if (conf >= FOREGROUND_CONFIDENCE) + knownFg.push_back(lab); + } + + /* + std::vector<CieLab> 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<unsigned int>(knownBg.size()), static_cast<unsigned int>(knownFg.size())); + + + std::vector<CieLab> 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<CieLab> 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<unsigned int, Tupel> hs; + + unsigned int progressResolution = pixelCount / 10; + + for (unsigned int i=0; i<pixelCount; i++) + { + if (i % progressResolution == 0) + { + float progress = + 30.0 + 60.0 * (float)i / (float)pixelCount; + //trace("### progress:%f", progress); + if (!progressReport(progress)) + { + error("User aborted"); + delete[] imageClab; + delete[] labelField; + workImage.setValid(false); + return workImage; + } + } + + if (cm[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<unsigned int, Tupel>::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<bgSignature.size() ; j++) + { + float d = lab.diffSq(bgSignature[j]); + if (d<minBg) + { + minBg = d; + minIndex = j; + } + } + Tupel tupel(0.0f, 0, 0.0f, 0); + tupel.minBgDist = minBg; + tupel.indexMinBg = minIndex; + float minFg = 1.0e6f; + minIndex = -1; + for (unsigned int j = 0 ; j < fgSignature.size() ; j++) + { + float d = lab.diffSq(fgSignature[j]); + if (d < minFg) + { + minFg = d; + minIndex = j; + } + } + tupel.minFgDist = minFg; + tupel.indexMinFg = minIndex; + if (fgSignature.empty()) + { + isBackground = (minBg <= clusterSize); + // remove next line to force behaviour of old algorithm + //error("foreground signature does not exist"); + //delete[] labelField; + //workImage.setValid(false); + //return workImage; + } + else + { + isBackground = minBg < minFg; + } + hs[image[i]] = tupel; + } + + if (isBackground) + cm[i] = CERTAIN_BACKGROUND_CONFIDENCE; + else + cm[i] = CERTAIN_FOREGROUND_CONFIDENCE; + } + } + + + delete[] imageClab; + + trace("### postProcessing"); + + + //#### postprocessing + smooth(cm, width, height, 0.333f, 0.333f, 0.333f); // average + normalizeMatrix(cm, pixelCount); + erode(cm, width, height); + keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.0/*sizeFactorToKeep*/); + + //for (int i=0; i < 2/*smoothness*/; i++) + // smooth(cm, width, height, 0.333f, 0.333f, 0.333f); // average + + normalizeMatrix(cm, pixelCount); + + for (unsigned int i=0; i < pixelCount; i++) + { + if (cm[i] >= 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<pixelCount ; i++) + { + float conf = cm[i]; + if (conf < FOREGROUND_CONFIDENCE) + image[i] = backgroundFillColor; + } + + delete[] labelField; + + trace("### Done"); + keepGoing = false; + return workImage; +} + + + +//############## +//## PRIVATE +//############## + +/** + * Initialize the Siox engine to its 'pristine' state. + * Performed at the beginning of extractForeground(). + */ +void Siox::init() +{ + limits[0] = 0.64f; + limits[1] = 1.28f; + limits[2] = 2.56f; + + float negLimits[3]; + negLimits[0] = -limits[0]; + negLimits[1] = -limits[1]; + negLimits[2] = -limits[2]; + + clusterSize = sqrEuclideanDist(limits, 3, negLimits); +} + + +/** + * Clean up any debris from processing. + */ +void Siox::cleanup() +{ +} + + + + +/** + * Stage 1 of the color signature work. 'dims' will be either + * 2 for grays, or 3 for colors + */ +void Siox::colorSignatureStage1(CieLab *points, + unsigned int leftBase, + unsigned int rightBase, + unsigned int recursionDepth, + unsigned int *clusterCount, + 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 + 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<CieLab> &inputVec, + std::vector<CieLab> &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<pixelCount ; idx++) + labelField[idx] = -1; + + int curlabel = 0; + int maxregion= 0; + int maxblob = 0; + + // slow but easy to understand: + std::vector<int> labelSizes; + for (unsigned long i=0 ; i<pixelCount ; i++) + { + int regionCount = 0; + if (labelField[i] == -1 && cm[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<pixelCount ; 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 Siox::depthFirstSearch(int startPos, + float threshold, int curLabel) +{ + // stores positions of labeled pixels, where the neighbours + // should still be checked for processing: + + //trace("startPos:%d threshold:%f curLabel:%d", + // startPos, threshold, curLabel); + + std::vector<int> 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<pixelCount ; idx++) + labelField[idx] = -1; + + //int maxRegion=0; // unused now + std::vector<int> pixelsToVisit; + for (unsigned long i=0; i<pixelCount; i++) + { // for all pixels + if (labelField[i]!=-1 || cm[i]<UNKNOWN_REGION_CONFIDENCE) + { + continue; // already visited or bg + } + + unsigned int origColor = image[i]; + unsigned long 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.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 + && 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; 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. + */ +void Siox::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]. + */ +void Siox::normalizeMatrix(float *cm, int cmSize) +{ + float max= -1000000.0f; + for (int i=0; i<cmSize; i++) + if (cm[i] > 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<cmSize; i++) + cm[i]=alpha*cm[i]; +} + +/** + * Blurs confidence matrix with a given symmetrically weighted kernel. + * + * In the standard case confidence matrix entries are between 0...1 and + * the weight factors sum up to 1. + */ +void Siox::smooth(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 Euclidean distance of p and q. + */ +float Siox::sqrEuclideanDist(float *p, int pSize, float *q) +{ + float sum=0.0; + for (int i=0; i<pSize; i++) + { + float v = p[i] - q[i]; + sum += v*v; + } + return sum; +} + + + + + + +} // namespace siox +} // namespace org + +//######################################################################## +//# E N D O F F I L E +//######################################################################## + |