PrusaSlicer-NonPlainar/xs/src/libslic3r/Point.cpp
bubnikv a6ea01a23f Moved some math macros (sqr, lerp, clamp) to libslic3r.h
Added UNUSED macro to libslic3r.h, used it to reduce some compile warnings.

Split the Int128 class from Clipper library to a separate file,
extended Int128 with intrinsic types wherever possible for performance,
added new geometric predicates.

Added a draft of new FillRectilinear3, which should reduce overfill near the perimeters in the future.
2017-07-27 10:39:43 +02:00

379 lines
9.4 KiB
C++

#include "Point.hpp"
#include "Line.hpp"
#include "MultiPoint.hpp"
#include <algorithm>
#include <cmath>
namespace Slic3r {
Point::Point(double x, double y)
{
this->x = lrint(x);
this->y = lrint(y);
}
std::string
Point::wkt() const
{
std::ostringstream ss;
ss << "POINT(" << this->x << " " << this->y << ")";
return ss.str();
}
std::string
Point::dump_perl() const
{
std::ostringstream ss;
ss << "[" << this->x << "," << this->y << "]";
return ss.str();
}
void
Point::scale(double factor)
{
this->x *= factor;
this->y *= factor;
}
void
Point::translate(double x, double y)
{
this->x += x;
this->y += y;
}
void
Point::translate(const Vector &vector)
{
this->translate(vector.x, vector.y);
}
void
Point::rotate(double angle)
{
double cur_x = (double)this->x;
double cur_y = (double)this->y;
double s = sin(angle);
double c = cos(angle);
this->x = (coord_t)round(c * cur_x - s * cur_y);
this->y = (coord_t)round(c * cur_y + s * cur_x);
}
void
Point::rotate(double angle, const Point &center)
{
double cur_x = (double)this->x;
double cur_y = (double)this->y;
double s = sin(angle);
double c = cos(angle);
double dx = cur_x - (double)center.x;
double dy = cur_y - (double)center.y;
this->x = (coord_t)round( (double)center.x + c * dx - s * dy );
this->y = (coord_t)round( (double)center.y + c * dy + s * dx );
}
bool
Point::coincides_with_epsilon(const Point &point) const
{
return std::abs(this->x - point.x) < SCALED_EPSILON && std::abs(this->y - point.y) < SCALED_EPSILON;
}
int
Point::nearest_point_index(const Points &points) const
{
PointConstPtrs p;
p.reserve(points.size());
for (Points::const_iterator it = points.begin(); it != points.end(); ++it)
p.push_back(&*it);
return this->nearest_point_index(p);
}
int Point::nearest_point_index(const PointConstPtrs &points) const
{
int idx = -1;
double distance = -1; // double because long is limited to 2147483647 on some platforms and it's not enough
for (PointConstPtrs::const_iterator it = points.begin(); it != points.end(); ++it) {
/* If the X distance of the candidate is > than the total distance of the
best previous candidate, we know we don't want it */
double d = sqr<double>(this->x - (*it)->x);
if (distance != -1 && d > distance) continue;
/* If the Y distance of the candidate is > than the total distance of the
best previous candidate, we know we don't want it */
d += sqr<double>(this->y - (*it)->y);
if (distance != -1 && d > distance) continue;
idx = it - points.begin();
distance = d;
if (distance < EPSILON) break;
}
return idx;
}
int
Point::nearest_point_index(const PointPtrs &points) const
{
PointConstPtrs p;
p.reserve(points.size());
for (PointPtrs::const_iterator it = points.begin(); it != points.end(); ++it)
p.push_back(*it);
return this->nearest_point_index(p);
}
bool
Point::nearest_point(const Points &points, Point* point) const
{
int idx = this->nearest_point_index(points);
if (idx == -1) return false;
*point = points.at(idx);
return true;
}
/* distance to the closest point of line */
double
Point::distance_to(const Line &line) const
{
const double dx = line.b.x - line.a.x;
const double dy = line.b.y - line.a.y;
const double l2 = dx*dx + dy*dy; // avoid a sqrt
if (l2 == 0.0) return this->distance_to(line.a); // line.a == line.b case
// Consider the line extending the segment, parameterized as line.a + t (line.b - line.a).
// We find projection of this point onto the line.
// It falls where t = [(this-line.a) . (line.b-line.a)] / |line.b-line.a|^2
const double t = ((this->x - line.a.x) * dx + (this->y - line.a.y) * dy) / l2;
if (t < 0.0) return this->distance_to(line.a); // beyond the 'a' end of the segment
else if (t > 1.0) return this->distance_to(line.b); // beyond the 'b' end of the segment
Point projection(
line.a.x + t * dx,
line.a.y + t * dy
);
return this->distance_to(projection);
}
double
Point::perp_distance_to(const Line &line) const
{
if (line.a.coincides_with(line.b)) return this->distance_to(line.a);
double n = (double)(line.b.x - line.a.x) * (double)(line.a.y - this->y)
- (double)(line.a.x - this->x) * (double)(line.b.y - line.a.y);
return std::abs(n) / line.length();
}
/* Three points are a counter-clockwise turn if ccw > 0, clockwise if
* ccw < 0, and collinear if ccw = 0 because ccw is a determinant that
* gives the signed area of the triangle formed by p1, p2 and this point.
* In other words it is the 2D cross product of p1-p2 and p1-this, i.e.
* z-component of their 3D cross product.
* We return double because it must be big enough to hold 2*max(|coordinate|)^2
*/
double
Point::ccw(const Point &p1, const Point &p2) const
{
return (double)(p2.x - p1.x)*(double)(this->y - p1.y) - (double)(p2.y - p1.y)*(double)(this->x - p1.x);
}
double
Point::ccw(const Line &line) const
{
return this->ccw(line.a, line.b);
}
// returns the CCW angle between this-p1 and this-p2
// i.e. this assumes a CCW rotation from p1 to p2 around this
double
Point::ccw_angle(const Point &p1, const Point &p2) const
{
double angle = atan2(p1.x - this->x, p1.y - this->y)
- atan2(p2.x - this->x, p2.y - this->y);
// we only want to return only positive angles
return angle <= 0 ? angle + 2*PI : angle;
}
Point
Point::projection_onto(const MultiPoint &poly) const
{
Point running_projection = poly.first_point();
double running_min = this->distance_to(running_projection);
Lines lines = poly.lines();
for (Lines::const_iterator line = lines.begin(); line != lines.end(); ++line) {
Point point_temp = this->projection_onto(*line);
if (this->distance_to(point_temp) < running_min) {
running_projection = point_temp;
running_min = this->distance_to(running_projection);
}
}
return running_projection;
}
Point
Point::projection_onto(const Line &line) const
{
if (line.a.coincides_with(line.b)) return line.a;
/*
(Ported from VisiLibity by Karl J. Obermeyer)
The projection of point_temp onto the line determined by
line_segment_temp can be represented as an affine combination
expressed in the form projection of
Point = theta*line_segment_temp.first + (1.0-theta)*line_segment_temp.second.
If theta is outside the interval [0,1], then one of the Line_Segment's endpoints
must be closest to calling Point.
*/
double lx = (double)(line.b.x - line.a.x);
double ly = (double)(line.b.y - line.a.y);
double theta = ( (double)(line.b.x - this->x)*lx + (double)(line.b.y- this->y)*ly )
/ ( sqr<double>(lx) + sqr<double>(ly) );
if (0.0 <= theta && theta <= 1.0)
return theta * line.a + (1.0-theta) * line.b;
// Else pick closest endpoint.
if (this->distance_to(line.a) < this->distance_to(line.b)) {
return line.a;
} else {
return line.b;
}
}
Point
Point::negative() const
{
return Point(-this->x, -this->y);
}
Vector
Point::vector_to(const Point &point) const
{
return Vector(point.x - this->x, point.y - this->y);
}
std::ostream&
operator<<(std::ostream &stm, const Pointf &pointf)
{
return stm << pointf.x << "," << pointf.y;
}
std::string
Pointf::wkt() const
{
std::ostringstream ss;
ss << "POINT(" << this->x << " " << this->y << ")";
return ss.str();
}
std::string
Pointf::dump_perl() const
{
std::ostringstream ss;
ss << "[" << this->x << "," << this->y << "]";
return ss.str();
}
void
Pointf::scale(double factor)
{
this->x *= factor;
this->y *= factor;
}
void
Pointf::translate(double x, double y)
{
this->x += x;
this->y += y;
}
void
Pointf::translate(const Vectorf &vector)
{
this->translate(vector.x, vector.y);
}
void
Pointf::rotate(double angle)
{
double cur_x = this->x;
double cur_y = this->y;
double s = sin(angle);
double c = cos(angle);
this->x = c * cur_x - s * cur_y;
this->y = c * cur_y + s * cur_x;
}
void
Pointf::rotate(double angle, const Pointf &center)
{
double cur_x = this->x;
double cur_y = this->y;
double s = sin(angle);
double c = cos(angle);
double dx = cur_x - center.x;
double dy = cur_y - center.y;
this->x = center.x + c * dx - s * dy;
this->y = center.y + c * dy + s * dx;
}
Pointf
Pointf::negative() const
{
return Pointf(-this->x, -this->y);
}
Vectorf
Pointf::vector_to(const Pointf &point) const
{
return Vectorf(point.x - this->x, point.y - this->y);
}
void
Pointf3::scale(double factor)
{
Pointf::scale(factor);
this->z *= factor;
}
void
Pointf3::translate(const Vectorf3 &vector)
{
this->translate(vector.x, vector.y, vector.z);
}
void
Pointf3::translate(double x, double y, double z)
{
Pointf::translate(x, y);
this->z += z;
}
double
Pointf3::distance_to(const Pointf3 &point) const
{
double dx = ((double)point.x - this->x);
double dy = ((double)point.y - this->y);
double dz = ((double)point.z - this->z);
return sqrt(dx*dx + dy*dy + dz*dz);
}
Pointf3
Pointf3::negative() const
{
return Pointf3(-this->x, -this->y, -this->z);
}
Vectorf3
Pointf3::vector_to(const Pointf3 &point) const
{
return Vectorf3(point.x - this->x, point.y - this->y, point.z - this->z);
}
}