https://github.com/davepagurek/StrokeStrip
Tip revision: c9ae2ebac9ecbc6461952df58a05288bccc571a1 authored by Dave Pagurek on 19 August 2021, 15:47:13 UTC
Add script to install and run on Mac
Add script to install and run on Mac
Tip revision: c9ae2eb
Polyline2D.h
#pragma once
#define _USE_MATH_DEFINES
#include <cmath>
#include <iostream>
#include <vector>
#include <cstdint>
#include <sstream>
#include <regex>
#include <assert.h>
#include <glm/glm.hpp>
extern double epsilon_small;
namespace SketchUI
{
inline std::vector<std::string> split_str(std::string str, std::string pattern) {
std::vector<std::string> strs;
size_t pos = 0;
std::string token;
while ((pos = str.find(pattern)) != std::string::npos) {
token = str.substr(0, pos);
strs.emplace_back(token);
str.erase(0, pos + pattern.length());
}
strs.emplace_back(str);
return strs;
}
class Point2D
{
public: // data
double x;
double y;
// the target tangent
// if none, these two value should be zeros
double tangent_x;
double tangent_y;
// curvature only used for preprocessing
double curvature;
public: // constructor
Point2D(void) : tangent_x(0), tangent_y(0), curvature(0) {}
Point2D(double p_x, double p_y) : x(p_x), y(p_y), tangent_x(0), tangent_y(0), curvature(0) {}
Point2D(double p_x, double p_y, double t_x, double t_y) : x(p_x), y(p_y), tangent_x(t_x), tangent_y(t_y), curvature(0) {}
Point2D(double p_x, double p_y, double t_x, double t_y, double curv) : x(p_x), y(p_y), tangent_x(t_x), tangent_y(t_y), curvature(curv) {}
public: // operations
const Point2D operator-(const Point2D& right) const
{
return Point2D(x - right.x, y - right.y);
}
const Point2D operator+(const Point2D& right) const
{
return Point2D(x + right.x, y + right.y);
}
Point2D& operator+=(const Point2D& right)
{
x += right.x;
y += right.y;
return *this;
}
Point2D& operator/=(double divider)
{
x /= divider;
y /= divider;
return *this;
}
Point2D& operator=(const Point2D& right)
{
x = right.x;
y = right.y;
tangent_x = right.tangent_x;
tangent_y = right.tangent_y;
return *this;
}
double SqLength(void) const { return x*x+y*y; }
double Length(void) const { return std::sqrt(x*x+y*y); }
void Normalize(void)
{
*this /= Length();
}
std::string to_string() const {
std::ostringstream oss;
oss << "\t" << x << "\t" << y;
return oss.str();
}
bool from_string(std::string point_str) {
point_str = point_str.substr(point_str.find_first_not_of("\t"));
std::vector<std::string> strs = SketchUI::split_str(point_str, "\t");
if (strs.size() < 2) {
return false;
}
x = std::stod(strs[0]);
y = std::stod(strs[1]);
return true;
}
};
class Polyline2D
{
public: // data
int stroke_ind;
int group_ind;
std::vector<std::pair<Point2D, std::int64_t>> points;
double width;
double sampling_interval;
public:
Polyline2D() : sampling_interval(5), width(epsilon_small) { }
// attributes
// calculate total length of the polyline as the sum of each line segment's length
double totalLen(void) const
{
if(points.size() == 0)
{
puts("Can't calculate polyline total length! There's no points.");
return -1.0;
}
else if(points.size() == 1)
{
return 0.0;
}
else
{
double len = 0.0;
for(size_t i = 0; i < points.size() - 1; i++)
{
len += (points[i+1].first - points[i].first).Length();
}
return len;
}
}
void computeSegLen(std::vector<double>& segLen) const
{
if(points.size() < 2)
{
puts("Can't calculate segment lengths! There's no segments.");
return;
}
segLen.resize(points.size() - 1);
for(size_t i = 0; i < points.size() - 1; i++)
{
segLen[i] = (points[i].first - points[i+1].first).Length();
}
}
void computeSegLen2(std::vector<double>& segLen) const
{
if(points.empty())
{
puts("Can't calculate segment lengths! There's no segments.");
return;
}
segLen.resize(points.size());
segLen[0] = 0;
for(size_t i = 1; i < points.size(); i++)
{
segLen[i] = (points[i].first - points[i - 1].first).Length();
}
}
void reparameterize(double sample_step_length, size_t min_samples = 3) {
double total = totalLen();
size_t num_samples = std::max(min_samples, (size_t)ceil(total / sample_step_length));
sample_step_length = total / num_samples;
std::vector<double> segLen;
computeSegLen(segLen);
sampling_interval = sample_step_length;
if (points.empty()) return;
std::vector<std::pair<Point2D, std::int64_t>> reparam_points;
reparam_points.reserve(num_samples + 2);
reparam_points.emplace_back(points[0]);
double acc = 0.;
double local_t;
for (size_t i = 0; i < segLen.size(); i++) {
acc += segLen[i];
double new_acc = acc;
size_t j = 1;
// reparameterize the current step
while (acc >= j * sample_step_length) {
local_t = (segLen[i] - (acc - j * sample_step_length)) / segLen[i];
Point2D reparam_p((1 - local_t) * points[i].first.x + local_t * points[i + 1].first.x,
(1 - local_t) * points[i].first.y + local_t * points[i + 1].first.y);
reparam_p.tangent_x = (1 - local_t) * points[i].first.tangent_x + local_t * points[i + 1].first.tangent_x;
reparam_p.tangent_y = (1 - local_t) * points[i].first.tangent_y + local_t * points[i + 1].first.tangent_y;
double length = std::sqrt(reparam_p.tangent_x * reparam_p.tangent_x + reparam_p.tangent_y * reparam_p.tangent_y);
if (length > std::numeric_limits<double>::epsilon()) {
reparam_p.tangent_x /= length;
reparam_p.tangent_y /= length;
}
long double reparam_t = (1 - local_t) * points[i].second + local_t * points[i + 1].second;
reparam_points.emplace_back(reparam_p, (std::int64_t)reparam_t);
// left segment length
new_acc = (acc - j * sample_step_length);
j++;
}
acc = new_acc;
}
if ((points[points.size() - 1].first - reparam_points.back().first).Length() > 1.) {
reparam_points.emplace_back(points[points.size() - 1]);
}
else {
reparam_points.back() = points[points.size() - 1];
}
points = reparam_points;
}
// calculate normalized tangent directions at each point,
// NOTE: there mustn't be duplicated points!!!
void getTangentDirs(std::vector<Point2D>& tangentDirs) const
{
tangentDirs.resize(0);
tangentDirs.reserve(points.size());
if(points.size() < 2)
{
puts("Can't get polyline tangent vectors! Polyline points less than 2.");
return;
}
for(size_t i = 0; i < points.size() - 1; i++)
{
tangentDirs.push_back(points[i+1].first-points[i].first);
tangentDirs.back().Normalize();
}
tangentDirs.push_back(tangentDirs.back());
for(size_t i = points.size() - 2; i > 0; i--)
{
tangentDirs[i] = tangentDirs[i-1] + tangentDirs[i];
tangentDirs[i].Normalize();
}
}
void reindex() {
if (points.size() == 0) return;
std::vector<std::pair<Point2D, std::int64_t>> dedup_points = points;
points.clear();
points.emplace_back(dedup_points.front());
for (size_t i = 1; i < dedup_points.size(); i++) {
if ((dedup_points[i].first - points.back().first).Length() > std::numeric_limits<double>::epsilon())
points.emplace_back(dedup_points[i]);
}
}
double get_t(size_t ind) const {
assert(ind < points.size());
std::vector<double> seg;
computeSegLen2(seg);
double all_seg = 0;
double total = totalLen();
for (auto &s : seg) {
all_seg += s;
s = all_seg / total;
}
return seg[ind];
}
bool is_roughly_circular() const {
//std::cout << "\t In is_roughly_circular" << std::endl;
if (points.empty()) return false;
if ((points.front().first - points.back().first).Length() < std::numeric_limits<double>::epsilon())
return true;
auto to_glm_vec = [](Point2D p)->glm::dvec2 {
return glm::dvec2(p.x, p.y);
};
double turning = 0;
double min_angle_dist = 5 * M_PI;
glm::dvec2 near_round_p = to_glm_vec(points.at(0).first);
for (size_t i = 1; i + 1 < points.size(); i++) {
glm::dvec2 t0 = to_glm_vec(points.at(i).first - points.at(i - 1).first);
glm::dvec2 t1 = to_glm_vec(points.at(i + 1).first - points.at(i).first);
t0 = glm::normalize(t0);
t1 = glm::normalize(t1);
glm::dvec2 up0(-t0.y, t0.x);
double sign = (glm::dot(up0, t1) > 0) ? 1 : -1;
turning += sign * std::acos(std::max(-1., std::min(glm::dot(t0, t1), 1.)));
double angle_dist = std::abs(std::abs(turning) - 2 * M_PI);
if (angle_dist < min_angle_dist) {
near_round_p = to_glm_vec(points.at(i).first);
min_angle_dist = angle_dist;
}
}
bool angle_check = (std::abs(turning) > 250. / 180 * M_PI);
bool dist_check = glm::distance(near_round_p, to_glm_vec(points.at(0).first)) < 20;
//std::cout << "\t Turn: " << turning/M_PI * 180 << "; " << glm::distance(near_round_p, to_glm_vec(points.at(0).first)) << std::endl;
return angle_check && dist_check;
}
void laplacian_smooth(size_t center_ind, double t_range) {
auto to_glm_vec = [](Point2D p)->glm::dvec2 {
return glm::dvec2(p.x, p.y);
};
std::vector<double> seg;
computeSegLen2(seg);
double all_seg = 0;
double total = totalLen();
for (auto &s : seg) {
all_seg += s;
s = all_seg / total;
}
std::vector<std::pair<Point2D, std::int64_t>> smooth_points = points;
for (int i = 0; i < smooth_points.size(); i++) {
if ((std::abs(seg[i] - seg[center_ind]) < t_range || std::abs(i - (int)center_ind) == 1)
&& i != 0 && i + 1 != smooth_points.size()) {
double wij = 1. / (points[i].first - points[i - 1].first).Length();//1
double wik = 1. / (points[i].first - points[i + 1].first).Length();//1
glm::dvec2 smoothed = (wij * to_glm_vec(points[i - 1].first) + wik * to_glm_vec(points[i + 1].first))
/ (wij + wik);
smooth_points[i].first = Point2D(smoothed.x, smoothed.y);
}
}
points = smooth_points;
}
void laplacian_smooth(size_t center_ind, size_t ind_range) {
auto to_glm_vec = [](Point2D p)->glm::dvec2 {
return glm::dvec2(p.x, p.y);
};
std::vector<std::pair<Point2D, std::int64_t>> smooth_points = points;
for (int i = 0; i < smooth_points.size(); i++) {
if ((std::abs(i - (int)center_ind) <= ind_range)
&& i != 0 && i + 1 != smooth_points.size()) {
double wij = 1. / (points[i].first - points[i - 1].first).Length();//1
double wik = 1. / (points[i].first - points[i + 1].first).Length();//1
glm::dvec2 smoothed = (wij * to_glm_vec(points[i - 1].first) + wik * to_glm_vec(points[i + 1].first))
/ (wij + wik);
smooth_points[i].first = Point2D(smoothed.x, smoothed.y);
}
}
points = smooth_points;
}
void laplacian_move(size_t center_ind, Point2D center_new, size_t ind_range, size_t num_itr = 1) {
auto to_glm_vec = [](Point2D p)->glm::dvec2 {
return glm::dvec2(p.x, p.y);
};
// First, compute laplacians before we move
struct Laplacian {
double lap_mag;
double angle;
Laplacian() : lap_mag(0), angle(0) { }
};
std::vector<Laplacian> laplacians;
laplacians.resize(points.size());
auto compute_laplacians = [&](std::vector<Laplacian> &laplacians) {
for (size_t i = 1; i + 1 < points.size(); i++) {
glm::dvec2 neighbor_vec = to_glm_vec(points[i + 1].first - points[i - 1].first);
glm::dvec2 neighbor_dir = glm::normalize(neighbor_vec);
glm::dvec2 lap = to_glm_vec(points[i].first) - 0.5 * to_glm_vec(points[i + 1].first + points[i - 1].first);
glm::dvec2 lap_dir = glm::normalize(lap);
glm::dvec2 neighbor_up(-neighbor_dir.y, neighbor_dir.x);
double sign = (glm::dot(neighbor_up, lap_dir) > 0) ? -1 : 1;
double dot_prod = std::max(-1., std::min(1., glm::dot(neighbor_dir, lap_dir)));
laplacians[i].angle = sign * std::acos(dot_prod);
laplacians[i].lap_mag = std::sqrt(glm::dot(lap, lap));
}
};
compute_laplacians(laplacians);
// Now we move
points[center_ind].first = center_new;
// Finally we smooth
for (size_t k = 0; k < num_itr; k++) {
std::vector<std::pair<Point2D, std::int64_t>> smooth_points = points;
for (int i = 0; i < smooth_points.size(); i++) {
if (i != 0 && i + 1 != smooth_points.size()) {
glm::dvec2 avg_neighbor = 0.5 * to_glm_vec(points[i + 1].first + points[i - 1].first);
glm::dvec2 neighbor_vec = to_glm_vec(points[i + 1].first - points[i - 1].first);
glm::dvec2 neighbor_dir = glm::normalize(neighbor_vec);
glm::dmat2x2 rotation(std::cos(laplacians[i].angle), -std::sin(laplacians[i].angle), std::sin(laplacians[i].angle), std::cos(laplacians[i].angle));
glm::dvec2 lap = rotation * neighbor_dir;
lap *= laplacians[i].lap_mag;
glm::dvec2 smoothed = avg_neighbor + lap;
if (!(std::abs(i - (int)center_ind) <= ind_range || std::abs(i - (int)center_ind) == 1)) {
smoothed = 0.1 * smoothed + 0.9 * to_glm_vec(points[i].first);
}
smooth_points[i].first = Point2D(smoothed.x, smoothed.y);
}
}
points = smooth_points;
}
}
Point2D get_endpoint_inward_tangent(size_t end_ind) const
{
Point2D tangent(0, 0);
if (end_ind == 0) {
tangent = points[1].first-points[0].first;
tangent.Normalize();
}
else if (end_ind == points.size() - 1) {
tangent = points[points.size() - 2].first-points[points.size() - 1].first;
tangent.Normalize();
}
return tangent;
}
// calculate center point of the polyline by equally weighted averaging
Point2D center(void) const
{
if(points.size() == 0)
{
puts("Can't calculate polyline center! There's no points.");
return Point2D(0.0, 0.0);
}
else if(points.size() == 1)
{
return points[0].first;
}
else
{
Point2D c(0.0, 0.0);
for(auto vi = points.begin(); vi != points.end(); vi++)
{
c += vi->first;
}
c.x /= points.size();
c.y /= points.size();
return c;
}
}
// compute normalized principal direction of the polyline by PCA
Point2D principalDir(void) const
{
if(points.size() < 2)
{
puts("Can't compute principal direction! Polyline points less than 2.");
return Point2D(0.0, 0.0);
}
else if(points.size() == 2)
{
return points[1].first - points[0].first;
}
else
{
Point2D centerPoint = this->center();
double a = 0.0, b = 0.0, d = 0.0;
for(size_t i = 0; i < points.size(); i++)
{
double x = points[i].first.x - centerPoint.x;
double y = points[i].first.y - centerPoint.y;
a += x * x;
b += x * y;
d += y * y;
}
double c = b;
double T = a + d, D = a*d - b*c;
double L1 = T/2 + sqrt(T*T/4-D);
if(c != 0.0)
{
Point2D dir(L1-d, c);
dir.Normalize();
return dir;
}
else if(b != 0.0)
{
Point2D dir(b, L1-a);
dir.Normalize();
return dir;
}
else
{
return Point2D(1.0, 0.0);
}
}
}
std::string to_string() const {
std::string poly_str = "{\n";
poly_str += "\t#" + std::to_string(stroke_ind) + "\t" + std::to_string(group_ind) + "\n";
poly_str += "\t@" + std::to_string(width) + "\n";
for (auto const & p : points) {
poly_str += p.first.to_string();
std::ostringstream oss;
oss << "\t" << p.second;
poly_str += oss.str() + "\n";
}
poly_str += "}\n";
return poly_str;
}
void from_string(std::string poly_str) {
stroke_ind = group_ind = -1;
poly_str = poly_str.substr(poly_str.find_first_not_of("{"));
std::vector<std::string> strs = SketchUI::split_str(poly_str, "\n");
width = epsilon_small;
for (size_t i = 0; i < strs.size(); i++) {
if (strs[i].empty())
continue;
if (strs[i].find("#") != std::string::npos) {
strs[i] = strs[i].substr(strs[i].find_first_of("#"));
assert(strs[i].length() > 1);
strs[i] = strs[i].substr(1);
std::vector<std::string> strs_digit = SketchUI::split_str(strs[i], "\t");
stroke_ind = std::stoi(strs_digit[0]);
group_ind = std::stoi(strs_digit[1]);
continue;
}
if (strs[i].find("@") != std::string::npos) {
strs[i] = strs[i].substr(strs[i].find_first_of("@"));
assert(strs[i].length() > 1);
strs[i] = strs[i].substr(1);
width = std::stoi(strs[i]);
continue;
}
Point2D point;
bool ok = point.from_string(strs[i]);
if (!ok) continue;
uint64_t time;
{
strs[i] = strs[i].substr(strs[i].find_first_not_of("\t"));
std::vector<std::string> strs_digit = SketchUI::split_str(strs[i], "\t");
time = std::stoul(strs_digit[2]);
}
// Deduplicate
if (points.size() == 0 ||
(points.back().first - point).SqLength() > std::numeric_limits<float>::epsilon()) {
Point2D prev;
if (points.size() > 0)
prev = points.back().first;
points.emplace_back(std::make_pair(point, time));
Point2D cur = points.back().first;
}
}
sampling_interval = 0;
for (size_t i = 0; i + 1 < points.size(); i++) {
sampling_interval += (points[i + 1].first - points[i].first).Length();
}
sampling_interval = (points.size() >= 2) ? sampling_interval/points.size() : 0;
}
std::string to_svg_string() const {
std::string svg_str = "";
if (points.empty()) return svg_str;
svg_str += "\t<path d=\"";
for (size_t i = 0; i < points.size(); i++) {
svg_str += (i == 0) ? "M " : "L ";
svg_str += std::to_string(points[i].first.x) + " " + std::to_string(points[i].first.y);
if (i + 1 < points.size())
svg_str += " ";
}
svg_str += "\" fill=\"none\" stroke=\"black\" stroke-width=\"" + std::to_string(this->width) + "\" />\n";
return svg_str;
}
};
}