jaseg/gerbolyze
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
gerbolyze/svg-flatten/src/nopencv.cpp
jaseg 4d711a2fe8 Update all svg-flatten dependencies
2023-10-26 00:32:02 +02:00

616 lines
20 KiB
C++
Raw Blame History

#include <iostream>
#include <iomanip>
#include <stack>
#include "nopencv.hpp"
#define STB_IMAGE_IMPLEMENTATION
#include <stb_image.h>
#define STB_IMAGE_RESIZE_IMPLEMENTATION
#include <stb_image_resize2.h>
#define IIR_GAUSS_BLUR_IMPLEMENTATION
#include "iir_gauss_blur.h"
template void iir_gauss_blur<uint8_t>(unsigned int width, unsigned int height, unsigned char components, uint8_t* image, float sigma);
template void iir_gauss_blur<uint32_t> (unsigned int width, unsigned int height, unsigned char components, uint32_t* image, float sigma);
template void iir_gauss_blur<float> (unsigned int width, unsigned int height, unsigned char components, float* image, float sigma);
using namespace gerbolyze;
using namespace gerbolyze::nopencv;
static constexpr bool debug = false;
/* directions:
* 0
* 7 1
* ^
* |
* 6 <--- X ---> 2
* |
* v
* 5 3
* 4
*
*/
enum Direction {
D_N,
D_NE,
D_E,
D_SE,
D_S,
D_SW,
D_W,
D_NW
};
const char * const dir_str[8] = { "N", "NE", "E", "SE", "S", "SW", "W", "NW" };
static struct {
int x;
int y;
} dir_to_coords[8] = {{0, -1}, {1, -1}, {1, 0}, {1, 1}, {0, 1}, {-1, 1}, {-1, 0}, {-1, -1}};
static Direction flip_direction[8] = {
D_S, /* 0 */
D_SW, /* 1 */
D_W, /* 2 */
D_NW, /* 3 */
D_N, /* 4 */
D_NE, /* 5 */
D_E, /* 6 */
D_SE /* 7 */
};
static void follow(gerbolyze::nopencv::Image32 &img, int start_x, int start_y, Direction initial_direction, int nbd, int connectivity, Polygon_i &poly) {
if (debug) {
cerr << "follow " << start_x << " " << start_y << " | dir=" << dir_str[initial_direction] << " nbd=" << nbd << " conn=" << connectivity << endl;
}
int dir_inc = (connectivity == 4) ? 2 : 1;
int probe_x, probe_y;
/* homing run: find starting point for algorithm steps below. */
bool found = false;
int k;
for (k=initial_direction; k<initial_direction+8; k += dir_inc) {
probe_x = start_x + dir_to_coords[k % 8].x;
probe_y = start_y + dir_to_coords[k % 8].y;
if (img.at_default(probe_x, probe_y) != 0) {
found = true;
break;
}
}
if (!found) { /* No nonzero pixels found. This is a single-pixel contour */
img.at(start_x, start_y) = nbd;
/* We must return these vertices counter-clockwise! */
poly.emplace_back(i2p{start_x, start_y+1});
poly.emplace_back(i2p{start_x+1, start_y+1});
poly.emplace_back(i2p{start_x+1, start_y});
poly.emplace_back(i2p{start_x, start_y});
return;
}
/* starting point found. */
int current_direction = k % 8;
int start_direction = current_direction;
int center_x = start_x, center_y = start_y;
if (debug) {
cerr << " init: " << center_x << " " << center_y << " / " << dir_str[current_direction] << endl;
}
do {
bool flag = false;
for (k = current_direction + 8 - dir_inc; k >= current_direction; k -= dir_inc) {
probe_x = center_x + dir_to_coords[k % 8].x;
probe_y = center_y + dir_to_coords[k % 8].y;
if (k%8 == D_E)
flag = true;
if (img.at_default(probe_x, probe_y) != 0) {
break;
}
}
int set_val = 0;
if (flag && img.at_default(center_x+1, center_y) == 0) {
img.at(center_x, center_y) = -nbd;
set_val = -nbd;
} else if (img.at(center_x, center_y) == 1) {
img.at(center_x, center_y) = nbd;
set_val = nbd;
}
for (int l = (current_direction + 8 - 2 + 1) / 2 * 2; l > k; l -= dir_inc) {
switch (l%8) {
case 0: poly.emplace_back(i2p{center_x, center_y}); break;
case 2: poly.emplace_back(i2p{center_x+1, center_y}); break;
case 4: poly.emplace_back(i2p{center_x+1, center_y+1}); break;
case 6: poly.emplace_back(i2p{center_x, center_y+1}); break;
}
}
center_x = probe_x;
center_y = probe_y;
current_direction = flip_direction[k % 8];
if (debug) {
cerr << " " << center_x << " " << center_y << " / " << dir_str[current_direction] << " -> " << set_val << endl;
}
} while (center_x != start_x || center_y != start_y || current_direction != start_direction);
}
void gerbolyze::nopencv::find_contours(gerbolyze::nopencv::Image32 &img, gerbolyze::nopencv::ContourCallback cb) {
/* Implementation of the hierarchical contour finding algorithm from Suzuki and Abe, 1983: Topological Structural
* Analysis of Digitized Binary Images by Border Following
*
* Written with these two resources as reference:
* https://theailearner.com/tag/suzuki-contour-algorithm-opencv/
* https://github.com/FreshJesh5/Suzuki-Algorithm/blob/master/contoursv1/contoursv1.cpp
*
* WARNING: input image MUST BE BINARIZE: All pixels must have value either 0 or 1. Otherwise, chaos ensues.
*/
int nbd = 1;
Polygon_i poly;
for (int y=0; y<img.rows(); y++) {
for (int x=0; x<img.cols(); x++) {
int val_xy = img.at(x, y);
/* Note: outer borders are followed with 8-connectivity, hole borders with 4-connectivity. This prevents
* incorrect results in this case:
*
* 1 1 1 | 0 0 0
* |
* 1 1 1 | 0 0 0
* ----------+---------- <== Here
* 0 0 0 | 1 1 1
* |
* 0 0 0 | 1 1 1
*/
if (img.at_default(x-1, y) == 0 && val_xy == 1) { /* outer border starting point */
nbd += 1;
follow(img, x, y, D_W, nbd, 8, poly);
cb(poly, CP_CONTOUR);
poly.clear();
} else if (val_xy >= 1 && img.at_default(x+1, y) == 0) { /* hole border starting point */
nbd += 1;
follow(img, x, y, D_E, nbd, 8, poly); /* FIXME should be 4? */
cb(poly, CP_HOLE);
poly.clear();
}
}
}
}
static size_t region_of_support(Polygon_i poly, size_t i) {
double x0 = poly[i][0], y0 = poly[i][1];
size_t sz = poly.size();
double last_l = 0;
double last_r = 0;
size_t k;
//cerr << "d: ";
for (k=1; k<(sz+1)/2; k++) {
size_t idx1 = (i + k) % sz;
size_t idx2 = (i + sz - k) % sz;
double x1 = poly[idx1][0], y1 = poly[idx1][1], x2 = poly[idx2][0], y2 = poly[idx2][1];
double l = sqrt(pow(x2-x1, 2) + pow(y2-y1, 2));
/* https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
* TODO: Check whether distance-to-line is an ok implementation here, the paper asks for distance to chord.
*/
double d = ((x2-x1)*(y1-y0) - (x1-x0)*(y2-y1)) / sqrt(pow(x2-x1, 2) + pow(y2-y1, 2));
//cerr << d << " ";
double r = d/l;
bool cond_a = l < last_l;
bool cond_b = ((d > 0) && (r < last_r)) || ((d < 0) && (r > last_r));
if (k > 2 && (cond_a || cond_b))
break;
last_l = l;
last_r = r;
}
//cerr << endl;
k -= 1;
return k;
}
int freeman_angle(const Polygon_i &poly, size_t i) {
/* f:
* 2
* 3 1
* ^
* |
* 4 <--- X ---> 0
* |
* v
* 5 7
* 6
*
*/
size_t sz = poly.size();
auto &p_last = poly[(i + sz - 1) % sz];
auto &p_now = poly[i];
auto dx = p_now[0] - p_last[0];
auto dy = p_now[1] - p_last[1];
/* both points must be neighbors */
assert (-1 <= dx && dx <= 1);
assert (-1 <= dy && dy <= 1);
assert (!(dx == 0 && dy == 0));
int lut[3][3] = {{3, 2, 1}, {4, -1, 0}, {5, 6, 7}};
return lut[dy+1][dx+1];
}
double k_curvature(const Polygon_i &poly, size_t i, size_t k) {
size_t sz = poly.size();
double acc = 0;
for (size_t idx = 0; idx < k; idx++) {
acc += freeman_angle(poly, (i + 2*sz - idx) % sz) - freeman_angle(poly, (i+idx + 1) % sz);
}
return acc / k;
}
double k_cos(const Polygon_i &poly, size_t i, size_t k) {
size_t sz = poly.size();
int64_t x0 = poly[i][0], y0 = poly[i][1];
int64_t x1 = poly[(i + sz + k) % sz][0], y1 = poly[(i + sz + k) % sz][1];
int64_t x2 = poly[(i + sz - k) % sz][0], y2 = poly[(i + sz - k) % sz][1];
auto xa = x0 - x1, ya = y0 - y1;
auto xb = x0 - x2, yb = y0 - y2;
auto dp = xa*yb + ya*xb;
auto sq_a = xa*xa + ya*ya;
auto sq_b = xb*xb + yb*yb;
return dp / (sqrt(sq_a)*sqrt(sq_b));
}
ContourCallback gerbolyze::nopencv::simplify_contours_teh_chin(ContourCallback cb) {
return [&cb](Polygon_i &poly, ContourPolarity cpol) {
size_t sz = poly.size();
vector<size_t> ros(sz);
vector<double> sig(sz);
vector<double> cur(sz);
vector<bool> retain(sz);
for (size_t i=0; i<sz; i++) {
ros[i] = region_of_support(poly, i);
sig[i] = fabs(k_cos(poly, i, ros[i]));
cur[i] = k_curvature(poly, i, 1);
retain[i] = true;
}
if (debug) {
cerr << endl;
cerr << "Polarity: " << cpol <<endl;
cerr << "Coords:"<<endl;
cerr << " x: ";
for (size_t i=0; i<sz; i++) {
cerr << setfill(' ') << setw(2) << poly[i][0] << " ";
}
cerr << endl;
cerr << " y: ";
for (size_t i=0; i<sz; i++) {
cerr << setfill(' ') << setw(2) << poly[i][1] << " ";
}
cerr << endl;
cerr << "Metrics:"<<endl;
cerr << "ros: ";
for (size_t i=0; i<sz; i++) {
cerr << setfill(' ') << setw(2) << ros[i] << " ";
}
cerr << endl;
cerr << "sig: ";
for (size_t i=0; i<sz; i++) {
cerr << setfill(' ') << setw(2) << sig[i] << " ";
}
cerr << endl;
}
/* Pass 0 (like opencv): Remove points with zero 1-curvature */
for (size_t i=0; i<sz; i++) {
if (cur[i] == 0) {
retain[i] = false;
break;
}
}
if (debug) {
cerr << "pass 0: ";
for (size_t i=0; i<sz; i++) {
cerr << (retain[i] ? "#" : ".");
}
cerr << endl;
}
/* 3a, Pass 1: Non-maxima suppression */
for (size_t i=0; i<sz; i++) {
for (size_t j=1; j<ros[i]/2; j++) {
if (sig[i] < sig[(i + j) % sz] || sig[i] < sig[(i + sz - j) % sz]) {
retain[i] = false;
break;
}
}
}
if (debug) {
cerr << "pass 1: ";
for (size_t i=0; i<sz; i++) {
cerr << (retain[i] ? "#" : ".");
}
cerr << endl;
}
/* 3b, Pass 2: Zero-curvature suppression */
for (size_t i=0; i<sz; i++) {
if (retain[i] && ros[i] == 1) {
if (sig[i] <= sig[(i + 1) % sz] || sig[i] <= sig[(i + sz - 1) % sz]) {
retain[i] = false;
}
}
}
if (debug) {
cerr << "pass 2: ";
for (size_t i=0; i<sz; i++) {
cerr << (retain[i] ? "#" : ".");
}
cerr << endl;
}
/* 3c, Pass 3: Further thinning */
for (size_t i=0; i<sz; i++) {
if (retain[i]) {
if (ros[i] == 1) {
if (retain[(i + sz - 1) % sz] || retain[(i + 1)%sz]) {
if (sig[i] < sig[(i + sz - 1)%sz] || sig[i] < sig[(i + 1)%sz]) {
retain[i] = false;
}
}
}
}
}
if (debug) {
cerr << "pass 3: ";
for (size_t i=0; i<sz; i++) {
cerr << (retain[i] ? "#" : ".");
}
cerr << endl;
}
Polygon_i new_poly;
for (size_t i=0; i<sz; i++) {
if (retain[i]) {
new_poly.push_back(poly[i]);
}
}
if (!new_poly.empty()) {
cb(new_poly, cpol);
}
};
}
static double dp_eps(double dx, double dy) {
/* Implementation of:
*
* Prasad, Dilip K., et al. "A novel framework for making dominant point detection methods non-parametric."
* Image and Vision Computing 30.11 (2012): 843-859.
* https://core.ac.uk/download/pdf/131287229.pdf
*
* For another implementation, see:
* https://github.com/BobLd/RamerDouglasPeuckerNetV2/blob/master/RamerDouglasPeuckerNetV2.Test/RamerDouglasPeuckerNetV2/RamerDouglasPeucker.cs
*/
double m = dy / dx;
double s = sqrt(pow(dx, 2) + pow(dy, 2));
double phi = atan(m);
double t_max = 1/s * (fabs(cos(phi)) + fabs(sin(phi)));
double t_max_polynomial = 1 - t_max + pow(t_max, 2);
return s * fmax(
atan(1/s * fabs(sin(phi) + cos(phi)) * t_max_polynomial),
atan(1/s * fabs(sin(phi) - cos(phi)) * t_max_polynomial));
}
/* a, b inclusive */
static array<size_t, 3> dp_step(Polygon_i &poly, size_t a, size_t b) {
double dx = poly[b][0] - poly[a][0];
double dy = poly[b][1] - poly[a][1];
double eps = dp_eps(dx, dy);
size_t max_idx = 0;
double max_dist = 0;
/* https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line */
double dist_ab = sqrt(pow(poly[b][0] - poly[a][0], 2) + pow(poly[b][1] - poly[a][1], 2));
for (size_t i=a+1; i<b; i++) {
double dist_i = fabs(
(poly[b][0] - poly[a][0]) * (poly[a][1] - poly[i][1])
- (poly[a][0] - poly[i][0]) * (poly[b][1] - poly[a][1]))
/ dist_ab;
if (dist_i > max_dist && dist_i > eps) {
max_dist = dist_i;
max_idx = i;
}
}
return {a, max_idx, b};
}
ContourCallback gerbolyze::nopencv::simplify_contours_douglas_peucker(ContourCallback cb) {
return [&cb](Polygon_i &poly, ContourPolarity cpol) {
Polygon_i out;
out.push_back(poly[0]);
stack<array<size_t, 3>> indices;
indices.push(dp_step(poly, 0, poly.size()-1));
while (!indices.empty()) {
auto idx = indices.top();
indices.pop(); /* awesome C++ api let's goooooo */
if (idx[1] > 0) {
indices.push(dp_step(poly, idx[0], idx[1]));
indices.push(dp_step(poly, idx[1], idx[2]));
} else {
out.push_back(poly[idx[2]]);
}
}
cb(out, cpol);
};
}
double gerbolyze::nopencv::polygon_area(Polygon_i &poly) {
double acc = 0;
size_t prev = poly.size() - 1;
for (size_t cur=0; cur<poly.size(); cur++) {
acc += (poly[prev][0] + poly[cur][0]) * (poly[prev][1] - poly[cur][1]);
prev = cur;
}
return acc / 2;
}
double gerbolyze::nopencv::polygon_perimeter(Polygon_i &poly) {
double acc = 0;
size_t prev = poly.size() - 1;
for (size_t cur=0; cur<poly.size(); cur++) {
double dx = poly[cur][0] - poly[prev][0];
double dy = poly[cur][1] - poly[prev][1];
acc += sqrt(dx*dx + dy*dy);
prev = cur;
}
return acc;
}
d2p gerbolyze::nopencv::polygon_centroid(Polygon_i &poly) {
double acc_x = 0, acc_y = 0;
double area = polygon_area(poly);
size_t prev = poly.size() - 1;
for (size_t cur=0; cur<poly.size(); cur++) {
double a = poly[prev][1]*poly[cur][0] - poly[cur][1]*poly[prev][0];
acc_x += (poly[prev][0] + poly[cur][0]) * a;
acc_y += (poly[prev][1] + poly[cur][1]) * a;
prev = cur;
}
return { acc_x / (6*area), acc_y / (6*area) };
}
template<typename T>
gerbolyze::nopencv::Image<T>::Image(int size_x, int size_y, const T *data) {
assert(size_x > 0 && size_x < 100000);
assert(size_y > 0 && size_y < 100000);
m_data = new T[size_x * size_y] { 0 };
m_rows = size_y;
m_cols = size_x;
if (data != nullptr) {
memcpy(m_data, data, sizeof(T) * size_x * size_y);
}
}
template<typename T>
bool gerbolyze::nopencv::Image<T>::load(const char *filename) {
return stb_to_internal(stbi_load(filename, &m_cols, &m_rows, nullptr, 1));
}
template<typename T>
bool gerbolyze::nopencv::Image<T>::load_memory(const void *buf, size_t len) {
return stb_to_internal(stbi_load_from_memory(reinterpret_cast<const uint8_t *>(buf), len, &m_cols, &m_rows, nullptr, 1));
}
template<typename T>
void gerbolyze::nopencv::Image<T>::binarize(T threshold) {
assert(m_data != nullptr);
assert(m_rows > 0 && m_cols > 0);
for (int y=0; y<m_rows; y++) {
for (int x=0; x<m_cols; x++) {
m_data[y*m_cols + x] = m_data[y*m_cols + x] >= threshold;
}
}
}
template<typename T>
bool gerbolyze::nopencv::Image<T>::stb_to_internal(uint8_t *data) {
if (data == nullptr)
return false;
if (m_rows < 0 || m_rows > 100000)
return false;
if (m_cols < 0 || m_cols > 100000)
return false;
m_data = new T[size()] { 0 };
for (int y=0; y<m_rows; y++) {
for (int x=0; x<m_cols; x++) {
m_data[y*m_cols + x] = data[y*m_cols + x];
}
}
stbi_image_free(data);
return true;
}
template<typename T>
void gerbolyze::nopencv::Image<T>::blur(int radius) {
iir_gauss_blur(m_cols, m_rows, 1, m_data, radius/2.0);
}
template<>
void gerbolyze::nopencv::Image<float>::resize(int new_w, int new_h) {
float *old_data = m_data;
m_data = new float[new_w * new_h];
stbir_resize_float_linear(old_data, m_cols, m_rows, 0,
m_data, new_w, new_h, 0,
STBIR_1CHANNEL);
m_cols = new_w;
m_rows = new_h;
delete old_data;
}
template<>
void gerbolyze::nopencv::Image<uint8_t>::resize(int new_w, int new_h) {
uint8_t *old_data = m_data;
m_data = new uint8_t[new_w * new_h];
stbir_resize_uint8_linear(old_data, m_cols, m_rows, 0,
m_data, new_w, new_h, 0,
STBIR_1CHANNEL);
m_cols = new_w;
m_rows = new_h;
delete old_data;
}
template gerbolyze::nopencv::Image<int32_t>::Image(int size_x, int size_y, const int32_t *data);
template bool gerbolyze::nopencv::Image<int32_t>::load(const char *filename);
template bool gerbolyze::nopencv::Image<int32_t>::load_memory(const void *buf, size_t len);
template void gerbolyze::nopencv::Image<int32_t>::binarize(int32_t threshold);
template bool gerbolyze::nopencv::Image<int32_t>::stb_to_internal(uint8_t *data);
template void gerbolyze::nopencv::Image<int32_t>::blur(int radius);
template gerbolyze::nopencv::Image<uint8_t>::Image(int size_x, int size_y, const uint8_t *data);
template bool gerbolyze::nopencv::Image<uint8_t>::load(const char *filename);
template bool gerbolyze::nopencv::Image<uint8_t>::load_memory(const void *buf, size_t len);
template void gerbolyze::nopencv::Image<uint8_t>::binarize(uint8_t threshold);
template bool gerbolyze::nopencv::Image<uint8_t>::stb_to_internal(uint8_t *data);
template void gerbolyze::nopencv::Image<uint8_t>::blur(int radius);
template gerbolyze::nopencv::Image<float>::Image(int size_x, int size_y, const float *data);
template bool gerbolyze::nopencv::Image<float>::load(const char *filename);
template bool gerbolyze::nopencv::Image<float>::load_memory(const void *buf, size_t len);
template void gerbolyze::nopencv::Image<float>::binarize(float threshold);
template bool gerbolyze::nopencv::Image<float>::stb_to_internal(uint8_t *data);
template void gerbolyze::nopencv::Image<float>::blur(int radius);
Reference in a new issue View git blame Copy permalink