![[tokkee]](http://tokkee.org/images/avatar.png){kind=avatar}
diff --git a/src/trace/siox.cpp b/src/trace/siox.cpp
index b0632567cfc65c4aaecdff4f8d2d31deee9a9e55..c69af04a165c85266a3897e98cdd8576018dac8e 100644 (file)
--- a/src/trace/siox.cpp
+++ b/src/trace/siox.cpp
#include <math.h>
#include <stdarg.h>
#include <map>
+#include <algorithm>
namespace org
//# C L A B
//########################################################################
-static std::map<unsigned long, CLAB> clabLookupTable;
-
/**
* Convert integer A, R, G, B values into an pixel value.
*/
-/**
- * Construct this CLAB from a packed-pixel ARGB value
- */
-CLAB::CLAB(unsigned long rgb)
+//#########################################
+//# 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)
{
- //First try looking up in the cache
- std::map<unsigned long, CLAB>::iterator iter;
- iter = clabLookupTable.find(rgb);
- if (iter != clabLookupTable.end())
+ double onetwo = x*qnrt(x);
+ return onetwo*onetwo;
+}
+
+
+static bool _clab_inited_ = false;
+void CieLab::init()
+{
+ if (!_clab_inited_)
{
- CLAB res = iter->second;
- C = res.C;
- L = res.L;
- A = res.A;
- B = res.B;
+ 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 a packed-pixel ARGB value
+ */
+CieLab::CieLab(unsigned long rgb)
+{
+ init();
+
int ir = (rgb>>16) & 0xff;
int ig = (rgb>> 8) & 0xff;
int ib = (rgb ) & 0xff;
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) 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) 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) pow((fb + 0.055) / 1.055, 2.4);
+ fb = (float) pow24((fb + 0.055) / 1.055);
else
fb = fb / 12.92;
- fr = fr * 100.0;
- fg = fg * 100.0;
- fb = fb * 100.0;
-
// Use white = D65
- float x = fr * 0.4124 + fg * 0.3576 + fb * 0.1805;
- float y = fr * 0.2126 + fg * 0.7152 + fb * 0.0722;
- float z = fr * 0.0193 + fg * 0.1192 + fb * 0.9505;
+ 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 / 95.047;
- float vy = y / 100.000;
- float vz = z / 108.883;
+ 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) 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) 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) pow(vz, 0.3333);
+ vz = (float) cbrt(vz);
else
vz = (7.787 * vz) + (16.0 / 116.0);
L = 116.0 * vy - 16.0;
A = 500.0 * (vx - vy);
B = 200.0 * (vy - vz);
-
- // Cache for next time
- clabLookupTable[rgb] = *this;
-
}
/**
- * Return this CLAB's value a a packed-pixel ARGB value
+ * Return this CieLab's value converted to a packed-pixel ARGB value
*/
-unsigned long CLAB::toRGB()
+unsigned long CieLab::toRGB()
{
float vy = (L + 16.0) / 116.0;
float vx = A / 500.0 + vy;
else
vz = (vz - 16.0 / 116.0) / 7.787;
- float x = 95.047 * vx; //use white = D65
- float y = 100.000 * vy;
- float z = 108.883 * vz;
-
- vx = x / 100.0;
- vy = y / 100.0;
- vz = z / 100.0;
+ 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);
/**
- * Computes squared euclidian distance in CLAB space for two colors
+ * Squared Euclidian 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 euclidian distance in CieLab space for two colors
* given as RGB values.
*/
-float CLAB::diffSq(unsigned int rgb1, unsigned int rgb2)
+float CieLab::diffSq(unsigned int rgb1, unsigned int rgb2)
{
- CLAB c1(rgb1);
- CLAB c2(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);
+ CieLab c1(rgb1);
+ CieLab c2(rgb2);
+ float euclid = c1.diffSq(c2);
return euclid;
}
/**
- * Computes squared euclidian distance in CLAB space for two colors
+ * Computes squared euclidian distance in CieLab space for two colors
* given as RGB values.
*/
-float CLAB::diff(unsigned int rgb0, unsigned int rgb1)
+float CieLab::diff(unsigned int rgb0, unsigned int rgb1)
{
return (float) sqrt(diffSq(rgb0, rgb1));
}
+
//########################################################################
//# T U P E L
//########################################################################
Tupel(float minBgDistArg, long indexMinBgArg,
float minFgDistArg, long indexMinFgArg)
{
- minBgDist = minBgDistArg;
- indexMinBg = indexMinBgArg;
- minFgDist = minFgDistArg;
- indexMinFg = 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;
+ 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;
+ minBgDist = other.minBgDist;
+ indexMinBg = other.indexMinBg;
+ minFgDist = other.minFgDist;
+ indexMinFg = other.indexMinFg;
+ return *this;
}
virtual ~Tupel()
{}
if (cmdata) delete[] cmdata;
}
+/**
+ * Error logging
+ */
+void SioxImage::error(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.
unsigned int y,
unsigned int pixval)
{
- if (x > width || y > height)
+ 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;
}
unsigned int g,
unsigned int b)
{
- if (x > width || y > height)
+ 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) |
*/
unsigned int SioxImage::getPixel(unsigned int x, unsigned int y)
{
- if (x > width || y > height)
+ 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];
}
unsigned int y,
float confval)
{
- if (x > width || y > height)
+ 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;
}
*/
float SioxImage::getConfidence(unsigned int x, unsigned int y)
{
- if (x > width || y > height)
+ if (x >= width || y >= height)
+ {
+ g_warning("getConfidence: out of bounds (%d,%d)/(%d,%d)",
+ x, y, width, height);
return 0.0;
+ }
unsigned long offset = width * y + x;
return cmdata[offset];
}
SioxImage Siox::extractForeground(const SioxImage &originalImage,
unsigned int backgroundFillColor)
{
+ trace("### Start");
+
init();
keepGoing = true;
trace("### Creating signatures");
//#### create color signatures
- std::vector<CLAB> knownBg;
- std::vector<CLAB> knownFg;
- std::vector<CLAB> imageClab;
+ 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);
- CLAB lab(pix);
+ CieLab lab(pix);
imageClab.push_back(lab);
if (cm <= BACKGROUND_CONFIDENCE)
- knownBg.push_back(lab); //note: uses CLAB(rgb)
+ 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:%d knownFg:%d", knownBg.size(), knownFg.size());
+ trace("knownBg:%zu knownFg:%zu", knownBg.size(), knownFg.size());
- std::vector<CLAB> bgSignature ;
+ std::vector<CieLab> bgSignature ;
if (!colorSignature(knownBg, bgSignature, 3))
{
error("Could not create background signature");
workImage.setValid(false);
+ delete[] imageClab;
delete[] labelField;
return workImage;
}
{
error("User aborted");
workImage.setValid(false);
+ delete[] imageClab;
delete[] labelField;
return workImage;
}
- std::vector<CLAB> fgSignature ;
+ std::vector<CieLab> fgSignature ;
if (!colorSignature(knownFg, fgSignature, 3))
{
error("Could not create foreground signature");
workImage.setValid(false);
+ delete[] imageClab;
delete[] labelField;
return workImage;
}
// segmentation impossible
error("Signature size is < 1. Segmentation is impossible");
workImage.setValid(false);
+ delete[] imageClab;
delete[] labelField;
return workImage;
}
{
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;
if (i % progressResolution == 0)
{
float progress =
- 30.0 + 60.0 * (float)i / (float)progressResolution;
+ 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;
std::map<unsigned int, Tupel>::iterator iter = hs.find(i);
if (iter != hs.end()) //found
{
- Tupel tupel = iter->second;
+ Tupel tupel = iter->second;
isBackground = tupel.minBgDist <= tupel.minFgDist;
}
else
{
- CLAB lab = imageClab[i];
- float minBg = sqrEuclidianDist(lab, bgSignature[0]);
- int minIndex=0;
+ CieLab lab = imageClab[i];
+ float minBg = lab.diffSq(bgSignature[0]);
+ int minIndex = 0;
for (unsigned int j=1; j<bgSignature.size() ; j++)
{
- float d = sqrEuclidianDist(lab, bgSignature[j]);
+ float d = lab.diffSq(bgSignature[j]);
if (d<minBg)
{
minBg = d;
Tupel tupel(0.0f, 0, 0.0f, 0);
tupel.minBgDist = minBg;
tupel.indexMinBg = minIndex;
- float minFg = 1.0e6f;
- minIndex = -1;
+ float minFg = 1.0e6f;
+ minIndex = -1;
for (unsigned int j = 0 ; j < fgSignature.size() ; j++)
{
- float d = sqrEuclidianDist(lab, fgSignature[j]);
+ float d = lab.diffSq(fgSignature[j]);
if (d < minFg)
{
minFg = d;
tupel.indexMinFg = minIndex;
if (fgSignature.size() == 0)
{
- isBackground=(minBg <= clusterSize);
+ isBackground = (minBg <= clusterSize);
// remove next line to force behaviour of old algorithm
//error("foreground signature does not exist");
//delete[] labelField;
}
+ delete[] imageClab;
trace("### postProcessing");
- //## postprocessing
- smooth(cm, width, height, 0.33f, 0.33f, 0.33f); // average
+ //#### 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.33f, 0.33f, 0.33f); // average
+ //for (int i=0; i < 2/*smoothness*/; i++)
+ // smooth(cm, width, height, 0.333f, 0.333f, 0.333f); // average
normalizeMatrix(cm, pixelCount);
}
- //#### Yaay. We are done. Now clear everything but the background
- for (unsigned int y = 0 ; y < height ; y++)
- for (unsigned int x = 0 ; x < width ; x++)
- {
- float conf = workImage.getConfidence(x, y);
- if (conf < FOREGROUND_CONFIDENCE)
- {
- workImage.setPixel(x, y, backgroundFillColor);
- }
- }
+ //#### 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;
* Stage 1 of the color signature work. 'dims' will be either
* 2 for grays, or 3 for colors
*/
-void Siox::colorSignatureStage1(CLAB *points,
+void Siox::colorSignatureStage1(CieLab *points,
unsigned int leftBase,
unsigned int rightBase,
unsigned int recursionDepth,
{
unsigned int currentDim = recursionDepth % dims;
- CLAB point = points[leftBase];
+ CieLab point = points[leftBase];
float min = point(currentDim);
float max = min;
else
{
//create a leaf
- CLAB newpoint;
+ CieLab newpoint;
newpoint.C = rightBase - leftBase;
/**
* Stage 2 of the color signature work
*/
-void Siox::colorSignatureStage2(CLAB *points,
+void Siox::colorSignatureStage2(CieLab *points,
unsigned int leftBase,
unsigned int rightBase,
unsigned int recursionDepth,
unsigned int currentDim = recursionDepth % dims;
- CLAB point = points[leftBase];
+ CieLab point = points[leftBase];
float min = point(currentDim);
float max = min;
if ((float)sum >= threshold)
{
float scale = (float)(rightBase - leftBase);
- CLAB newpoint;
+ CieLab newpoint;
for (; leftBase < rightBase; leftBase++)
newpoint.add(points[leftBase]);
/**
* Main color signature method
*/
-bool Siox::colorSignature(const std::vector<CLAB> &inputVec,
- std::vector<CLAB> &result,
+bool Siox::colorSignature(const std::vector<CieLab> &inputVec,
+ std::vector<CieLab> &result,
const unsigned int dims)
{
if (length < 1) // no error. just don't do anything
return true;
- CLAB *input = (CLAB *) malloc(length * sizeof(CLAB));
+ CieLab *input = (CieLab *) malloc(length * sizeof(CieLab));
if (!input)
{
// check all four neighbours
int left = pos-1;
if (((int)x)-1 >= 0 && labelField[left] == -1
- && CLAB::diff(image[left], origColor)<1.0)
+ && CieLab::diff(image[left], origColor)<1.0)
{
labelField[left]=curLabel;
cm[left]=CERTAIN_FOREGROUND_CONFIDENCE;
}
int right = pos+1;
if (x+1 < width && labelField[right]==-1
- && CLAB::diff(image[right], origColor)<1.0)
+ && CieLab::diff(image[right], origColor)<1.0)
{
labelField[right]=curLabel;
cm[right]=CERTAIN_FOREGROUND_CONFIDENCE;
}
int top = pos - width;
if (((int)y)-1>=0 && labelField[top]==-1
- && CLAB::diff(image[top], origColor)<1.0)
+ && CieLab::diff(image[top], origColor)<1.0)
{
labelField[top]=curLabel;
cm[top]=CERTAIN_FOREGROUND_CONFIDENCE;
}
int bottom = pos + width;
if (y+1 < height && labelField[bottom]==-1
- && CLAB::diff(image[bottom], origColor)<1.0)
+ && CieLab::diff(image[bottom], origColor)<1.0)
{
labelField[bottom]=curLabel;
cm[bottom]=CERTAIN_FOREGROUND_CONFIDENCE;
-
/**
* Normalizes the matrix to values to [0..1].
*/
{
float max= -1000000.0f;
for (int i=0; i<cmSize; i++)
- if (max<cm[i] > max)
- max=cm[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;
return sum;
}
-/**
- * Squared Euclidian distance of p and q.
- */
-float Siox::sqrEuclidianDist(const CLAB &p, const CLAB &q)
-{
- float sum=0;
- 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;
-}
-
-
-