From c853ffb5b2f75f5a889ed2e3ef89b818a736e87a Mon Sep 17 00:00:00 2001 From: Daniel Baumann Date: Sat, 13 Apr 2024 13:50:49 +0200 Subject: Adding upstream version 1.3+ds. Signed-off-by: Daniel Baumann --- src/trace/siox.cpp | 770 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 770 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..a958821 --- /dev/null +++ b/src/trace/siox.cpp @@ -0,0 +1,770 @@ +// 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 +#include +#include +#include +#include +#include + +#include "siox.h" +#include "async/progress.h" + +namespace Inkscape { +namespace Trace { + +//######################################################################## +//# S I O X I M A G E +//######################################################################## + +SioxImage::SioxImage(Glib::RefPtr const &buf) +{ + width = buf->get_width(); + height = buf->get_height(); + + // Allocate data arrays. + int size = width * height; + pixdata.resize(size); + cmdata.resize(size); + + int rowstride = buf->get_rowstride(); + int nchannels = buf->get_n_channels(); + auto data = buf->get_pixels(); + + // Copy pixel data. + for (int y = 0; y < height; y++) { + auto p = data + rowstride * y; + for (int x = 0; x < width; x++) { + uint32_t r = p[0]; + uint32_t g = p[1]; + uint32_t b = p[2]; + uint32_t a = nchannels == 3 ? 255 : p[3]; + pixdata[offset(x, y)] = (a << 24) | (r << 16) | (g << 8) | b; + p += nchannels; + } + } + + // Zero confidence matrix. + std::fill(cmdata.begin(), cmdata.end(), 0.0f); +} + +Glib::RefPtr SioxImage::getGdkPixbuf() const +{ + auto buf = Gdk::Pixbuf::create(Gdk::COLORSPACE_RGB, true, 8, width, height); + + int rowstride = buf->get_rowstride(); + int nchannels = buf->get_n_channels(); + auto data = buf->get_pixels(); + + for (int y = 0; y < height; y++) { + auto p = data + rowstride * y; + for (int x = 0; x < width; x++) { + uint32_t rgb = pixdata[offset(x, y)]; + p[0] = (rgb >> 16) & 0xff; // r + p[1] = (rgb >> 8) & 0xff; // g + p[2] = (rgb ) & 0xff; // b + p[3] = (rgb >> 24) & 0xff; // a + p += nchannels; + } + } + + return buf; +} + +bool SioxImage::writePPM(char const *filename) const +{ + auto f = std::fopen(filename, "wb"); + if (!f) { + return false; + } + + std::fprintf(f, "P6 %u %u 255\n", width, height); + + for (int y = 0; y < height; y++) { + for (int x = 0; x < width; x++) { + uint32_t rgb = pixdata[offset(x, y)]; + uint8_t r = (rgb >> 16) & 0xff; + uint8_t g = (rgb >> 8) & 0xff; + uint8_t b = (rgb ) & 0xff; + std::fputc(r, f); + std::fputc(g, f); + std::fputc(b, f); + } + } + + std::fclose(f); + + return true; +} + +unsigned SioxImage::hash() const +{ + unsigned result = width * height; + + for (int i = 0; i < width * height; i++) { + result = 3 * result + (unsigned)pixdata[i] + (unsigned)((1 << 16) * cmdata[i]); + } + + return result; +} + +//######################################################################## +//# S I O X +//######################################################################## + +namespace { + +/** + * Apply a function which updates each pixel depending on the value of its neighbours. + */ +template +void apply_adjacent(float *cm, int xres, int yres, F f) +{ + for (int y = 0; y < yres; y++) { + for (int x = 0; x < xres - 1; x++) { + int idx = y * xres + x; + f(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; + f(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; + f(cm[idx], cm[idx + xres]); + } + } + for (int y = yres - 1; y >= 1; y--) { + for (int x = 0; x < xres; x++) { + int idx = y * xres + x; + f(cm[idx], cm[idx - xres]); + } + } +} + +/** + * Applies the morphological dilate operator. + * + * Can be used to close small holes in the given confidence matrix. + */ +void dilate(float *cm, int xres, int yres) +{ + apply_adjacent(cm, xres, yres, [] (float &a, float b) { + if (b > a) { + a = b; + } + }); +} + +/** + * Applies the morphological erode operator. + */ +void erode(float *cm, int xres, int yres) +{ + apply_adjacent(cm, xres, yres, [] (float &a, float b) { + if (b < a) { + a = b; + } + }); +} + +/** + * Multiplies matrix with the given scalar. + */ +void premultiplyMatrix(float alpha, float *cm, int cmSize) +{ + for (int i = 0; i < cmSize; i++) { + cm[i] *= alpha; + } +} + +/** + * Normalizes the matrix to values to [0..1]. + */ +void normalizeMatrix(float *cm, int cmSize) +{ + float max = 0.0f; + for (int i = 0; i < cmSize; i++) { + if (cm[i] > max) { + max = cm[i]; + } + } + + if (max <= 0.0f || max == 1.0f) { + return; + } + + float alpha = 1.0f / max; + premultiplyMatrix(alpha, cm, cmSize); +} + +/** + * 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 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 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::Siox(Async::Progress &progress) + : progress(&progress) + , width(0) + , height(0) + , pixelCount(0) + , image(nullptr) + , cm(nullptr) {} + +void Siox::error(std::string const &msg) +{ + g_warning("Siox error: %s\n", msg.c_str()); +} + +void Siox::trace(std::string const &msg) +{ + g_message("Siox: %s\n", msg.c_str()); +} + +SioxImage Siox::extractForeground(SioxImage const &originalImage, uint32_t backgroundFillColor) +{ + trace("### Start"); + + init(); + + SioxImage workImage = originalImage; + + // Fetch some info from the image. + width = workImage.getWidth(); + height = workImage.getHeight(); + pixelCount = width * height; + image = workImage.getImageData(); + cm = workImage.getConfidenceData(); + + // Create labelField. + auto labelField_storage = std::make_unique(pixelCount); + labelField = labelField_storage.get(); + + trace("### Creating signatures"); + + // Create color signatures. + std::vector knownBg, knownFg; + auto imageClab = std::make_unique(pixelCount); + for (int i = 0; i < pixelCount; i++) { + float conf = cm[i]; + uint32_t pix = image[i]; + CieLab lab = pix; + imageClab[i] = lab; + if (conf <= BACKGROUND_CONFIDENCE) { + knownBg.emplace_back(lab); + } else if (conf >= FOREGROUND_CONFIDENCE) { + knownFg.emplace_back(lab); + } + } + + progress->report_or_throw(0.1); + + trace("knownBg:" + std::to_string(knownBg.size()) + " knownFg:" + std::to_string(knownFg.size())); + + std::vector bgSignature; + colorSignature(knownBg, bgSignature, 3); + + progress->report_or_throw(0.2); + + std::vector fgSignature; + colorSignature(knownFg, fgSignature, 3); + + // trace("### bgSignature:" + std::to_string(bgSignature.size())); + + if (bgSignature.empty()) { + // segmentation impossible + error("Signature size is < 1. Segmentation is impossible"); + throw Exception(); + } + + progress->report_or_throw(0.3); + + // classify using color signatures, + // classification cached in hashmap for drb and speedup purposes + trace("### Analyzing image"); + + std::unordered_map hs; + + int progressResolution = pixelCount / 10; + + for (int i = 0; i < pixelCount; i++) { + if (i % progressResolution == 0) { + progress->report_or_throw(0.3 + 0.6 * i / pixelCount); + } + + 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 + auto [it, inserted] = hs.emplace(image[i], false); + if (inserted) { + auto const &lab = imageClab[i]; + + float minBg = std::numeric_limits::max(); + for (auto const &s : bgSignature) { + minBg = std::min(minBg, CieLab::diffSq(lab, s)); + } + + float minFg; + if (fgSignature.empty()) { + minFg = clusterSize; + } else { + minFg = std::numeric_limits::max(); + for (auto const &s : fgSignature) { + minFg = std::min(minFg, CieLab::diffSq(lab, s)); + } + } + + it->second = minBg < minFg; + } + + bool isBackground = it->second; + cm[i] = isBackground ? CERTAIN_BACKGROUND_CONFIDENCE : CERTAIN_FOREGROUND_CONFIDENCE; + } + } + + hs.clear(); + imageClab.reset(); + + 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 (int i = 0; i < pixelCount; i++) { + cm[i] = cm[i] >= UNKNOWN_REGION_CONFIDENCE + ? CERTAIN_FOREGROUND_CONFIDENCE + : CERTAIN_BACKGROUND_CONFIDENCE; + } + + keepOnlyLargeComponents(UNKNOWN_REGION_CONFIDENCE, 1.5/*sizeFactorToKeep*/); + fillColorRegions(); + dilate(cm, width, height); + + progress->report_or_throw(1.0); + + // We are done. Now clear everything but the background. + for (int i = 0; i < pixelCount; i++) { + if (cm[i] < FOREGROUND_CONFIDENCE) { + image[i] = backgroundFillColor; + } + } + + trace("### Done"); + return workImage; +} + +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); +} + +void Siox::colorSignatureStage1(CieLab *points, + unsigned leftBase, + unsigned rightBase, + unsigned recursionDepth, + unsigned *clusterCount, + unsigned dims) +{ + unsigned currentDim = recursionDepth % dims; + CieLab point = points[leftBase]; + float min = point(currentDim); + float max = min; + + for (unsigned 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 left = leftBase; + unsigned 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.0f / newpoint.C); + } + points[*clusterCount] = newpoint; + (*clusterCount)++; + } +} + +void Siox::colorSignatureStage2(CieLab *points, + unsigned leftBase, + unsigned rightBase, + unsigned recursionDepth, + unsigned *clusterCount, + float threshold, + unsigned dims) +{ + unsigned currentDim = recursionDepth % dims; + CieLab point = points[leftBase]; + float min = point(currentDim); + float max = min; + + for (unsigned 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 left = leftBase; + unsigned 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 sum = 0; + for (unsigned i = leftBase; i < rightBase; i++) { + sum += points[i].C; + } + + if (sum >= threshold) { + float scale = 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)++; + } + } +} + +void Siox::colorSignature(std::vector const &inputVec, + std::vector &result, + unsigned dims) +{ + if (inputVec.empty()) { // no error. just don't do anything + return; + } + + unsigned length = inputVec.size(); + result = inputVec; + + unsigned stage1length = 0; + colorSignatureStage1(result.data(), 0, length, 0, &stage1length, dims); + + unsigned stage2length = 0; + colorSignatureStage2(result.data(), 0, stage1length, 0, &stage2length, length * 0.001, dims); + + result.resize(stage2length); +} + +void Siox::keepOnlyLargeComponents(float threshold, double sizeFactorToKeep) +{ + for (int idx = 0; idx < pixelCount; idx++) { + labelField[idx] = -1; + } + + int curlabel = 0; + int maxregion = 0; + int maxblob = 0; + + // slow but easy to understand: + std::vector labelSizes; + for (int i = 0; i < pixelCount; i++) { + int regionCount = 0; + if (labelField[i] == -1 && cm[i] >= threshold) { + regionCount = depthFirstSearch(i, threshold, curlabel++); + labelSizes.emplace_back(regionCount); + } + + if (regionCount > maxregion) { + maxregion = regionCount; + maxblob = curlabel-1; + } + } + + for (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 pixelsToVisit; + int componentSize = 0; + + if (labelField[startPos] == -1 && cm[startPos] >= threshold) { + labelField[startPos] = curLabel; + componentSize++; + pixelsToVisit.emplace_back(startPos); + } + + while (!pixelsToVisit.empty()) { + int pos = pixelsToVisit[pixelsToVisit.size() - 1]; + pixelsToVisit.erase(pixelsToVisit.end() - 1); + int x = pos % width; + int y = pos / width; + + // check all four neighbours + int left = pos - 1; + if (x - 1 >= 0 && labelField[left] == -1 && cm[left] >= threshold) { + labelField[left] = curLabel; + componentSize++; + pixelsToVisit.emplace_back(left); + } + + int right = pos + 1; + if (x + 1 < width && labelField[right] == -1 && cm[right] >= threshold) { + labelField[right] = curLabel; + componentSize++; + pixelsToVisit.emplace_back(right); + } + + int top = pos - width; + if (y - 1 >= 0 && labelField[top] == -1 && cm[top] >= threshold) { + labelField[top] = curLabel; + componentSize++; + pixelsToVisit.emplace_back(top); + } + + int bottom = pos + width; + if (y + 1 < height && labelField[bottom] == -1 && cm[bottom] >= threshold) { + labelField[bottom] = curLabel; + componentSize++; + pixelsToVisit.emplace_back(bottom); + } + } + + return componentSize; +} + +void Siox::fillColorRegions() +{ + for (int idx = 0; idx < pixelCount; idx++) { + labelField[idx] = -1; + } + + std::vector pixelsToVisit; + for (int i = 0; i < pixelCount; i++) { // for all pixels + if (labelField[i] != -1 || cm[i] < UNKNOWN_REGION_CONFIDENCE) { + continue; // already visited or bg + } + + uint32_t origColor = image[i]; + int curLabel = i+1; + labelField[i] = curLabel; + cm[i] = CERTAIN_FOREGROUND_CONFIDENCE; + + // int componentSize = 1; + pixelsToVisit.emplace_back(i); + // depth first search to fill region + while (!pixelsToVisit.empty()) { + int pos = pixelsToVisit[pixelsToVisit.size() - 1]; + pixelsToVisit.erase(pixelsToVisit.end() - 1); + int x = pos % width; + int y = pos / width; + // check all four neighbours + int left = pos - 1; + if (x - 1 >= 0 && labelField[left] == -1 && CieLab::diff(image[left], origColor) < 1.0) { + labelField[left] = curLabel; + cm[left] = CERTAIN_FOREGROUND_CONFIDENCE; + // ++componentSize; + pixelsToVisit.emplace_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.emplace_back(right); + } + int top = pos - width; + if (y - 1 >= 0 && labelField[top] == -1 && CieLab::diff(image[top], origColor) < 1.0) { + labelField[top] = curLabel; + cm[top] = CERTAIN_FOREGROUND_CONFIDENCE; + // ++componentSize; + pixelsToVisit.emplace_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.emplace_back(bottom); + } + } + } +} + +} // namespace Trace +} // namespace Inkscape -- cgit v1.2.3