981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <algorithm>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include "helper/geom.h"

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include "helper/geom-curves.h"

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <typeinfo>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/pathvector.h>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/path.h>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/bezier-curve.h>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/hvlinesegment.h>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/transforms.h>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/rect.h>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/coord.h>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <2geom/sbasis-to-bezier.h>

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen#include <math.h> // for M_PI

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelenusing Geom::X;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelenusing Geom::Y;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

40742313779ee5e43be93a9191f1c86412cf183bKrzysztof Kosiński

8001ba81cb851b38d86650a2fef5817facffb763johanengelenstatic void

8001ba81cb851b38d86650a2fef5817facffb763johanengelencubic_bbox (Geom::Coord x000, Geom::Coord y000, Geom::Coord x001, Geom::Coord y001, Geom::Coord x011, Geom::Coord y011, Geom::Coord x111, Geom::Coord y111, Geom::Rect &bbox)

8001ba81cb851b38d86650a2fef5817facffb763johanengelen{

40742313779ee5e43be93a9191f1c86412cf183bKrzysztof Kosiński Geom::Coord a, b, c, D;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński

29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen bbox[0].expandTo(x111);

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński bbox[1].expandTo(y111);

6c3e745a94ef6b25a4ef9f018d350a7535aa45afTed Gould

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński // It already contains (x000,y000) and (x111,y111)

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen // All points of the Bezier lie in the convex hull of (x000,y000), (x001,y001), (x011,y011) and (x111,y111)

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen // So, if it also contains (x001,y001) and (x011,y011) we don't have to compute anything else!

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen // Note that we compute it for the X and Y range separately to make it easier to use them below

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bool containsXrange = bbox[0].contains(x001) && bbox[0].contains(x011);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bool containsYrange = bbox[1].contains(y001) && bbox[1].contains(y011);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /*

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (!containsXrange) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen a = 3 * x000 - 9 * x001 + 9 * x011 - 3 * x111;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen b = 6 * x001 - 12 * x011 + 6 * x111;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen c = 3 * x011 - 3 * x111;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen /*

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (fabs (a) < Geom::EPSILON) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* s = -c / b */

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (fabs (b) > Geom::EPSILON) {

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double s;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen s = -c / b;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if ((s > 0.0) && (s < 1.0)) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen double t = 1.0 - s;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński bbox[0].expandTo(xttt);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen }

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen }

6492fcbc5fbe9a07962f989247731b6435fd72b9johanengelen } else {

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński D = b * b - 4 * a * c;

29684a16b6c92bee28a94fdc2607bcc143950fa8johanengelen if (D >= 0.0) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Coord d, s, t, xttt;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen /* Have solution */

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen d = sqrt (D);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen s = (-b + d) / (2 * a);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if ((s > 0.0) && (s < 1.0)) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen t = 1.0 - s;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bbox[0].expandTo(xttt);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen }

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen s = (-b - d) / (2 * a);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if ((s > 0.0) && (s < 1.0)) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen t = 1.0 - s;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen xttt = s * s * s * x000 + 3 * s * s * t * x001 + 3 * s * t * t * x011 + t * t * t * x111;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bbox[0].expandTo(xttt);

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński }

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński }

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński }

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński }

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński if (!containsYrange) {

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński a = 3 * y000 - 9 * y001 + 9 * y011 - 3 * y111;

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński b = 6 * y001 - 12 * y011 + 6 * y111;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński c = 3 * y011 - 3 * y111;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (fabs (a) < Geom::EPSILON) {

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński /* s = -c / b */

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński if (fabs (b) > Geom::EPSILON) {

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double s;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński s = -c / b;

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński if ((s > 0.0) && (s < 1.0)) {

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double t = 1.0 - s;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński double yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński bbox[1].expandTo(yttt);

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński }

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński }

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński } else {

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński /* s = (-b +/- sqrt (b * b - 4 * a * c)) / 2 * a; */

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński D = b * b - 4 * a * c;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if (D >= 0.0) {

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński Geom::Coord d, s, t, yttt;

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński /* Have solution */

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński d = sqrt (D);

e4369b05aaa20df73a37f4d319ce456865cc74f3Krzysztof Kosiński s = (-b + d) / (2 * a);

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński if ((s > 0.0) && (s < 1.0)) {

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński t = 1.0 - s;

76addc201c409e81eaaa73fe27cc0f79c4db097cKrzysztof Kosiński yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bbox[1].expandTo(yttt);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen }

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen s = (-b - d) / (2 * a);

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński if ((s > 0.0) && (s < 1.0)) {

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński t = 1.0 - s;

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński yttt = s * s * s * y000 + 3 * s * s * t * y001 + 3 * s * t * t * y011 + t * t * t * y111;

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński bbox[1].expandTo(yttt);

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński }

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński }

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński }

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński }

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński}

00f9ca0b3aa57e09f3c3f3632c5427fc03499df5Krzysztof Kosiński

981b809bc6ed10a21e89444d9447e5475801874fjohanengelenGeom::OptRect

981b809bc6ed10a21e89444d9447e5475801874fjohanengelenbounds_fast_transformed(Geom::PathVector const & pv, Geom::Affine const & t)

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen{

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return bounds_exact_transformed(pv, t); //use this as it is faster for now! :)

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen}

e6bdf746e2d9e775704a475a29cc1bb167ec271cjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelenGeom::OptRect

981b809bc6ed10a21e89444d9447e5475801874fjohanengelenbounds_exact_transformed(Geom::PathVector const & pv, Geom::Affine const & t)

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen{

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (pv.empty())

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen return Geom::OptRect();

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Point initial = pv.front().initialPoint() * t;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Rect bbox(initial, initial); // obtain well defined bbox as starting point to unionWith

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (Geom::PathVector::const_iterator it = pv.begin(); it != pv.end(); ++it) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bbox.expandTo(it->initialPoint() * t);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen // don't loop including closing segment, since that segment can never increase the bbox

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen for (Geom::Path::const_iterator cit = it->begin(); cit != it->end_open(); ++cit) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Curve const &c = *cit;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen unsigned order = 0;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (Geom::BezierCurve const* b = dynamic_cast<Geom::BezierCurve const*>(&c)) {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen order = b->order();

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen }

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen if (order == 1) { // line segment

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bbox.expandTo(c.finalPoint() * t);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

40742313779ee5e43be93a9191f1c86412cf183bKrzysztof Kosiński // TODO: we can make the case for quadratics faster by degree elevating them to

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen // cubic and then taking the bbox of that.

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen } else if (order == 3) { // cubic bezier

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::CubicBezier const &cubic_bezier = static_cast<Geom::CubicBezier const&>(c);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Point c0 = cubic_bezier[0] * t;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Point c1 = cubic_bezier[1] * t;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Point c2 = cubic_bezier[2] * t;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Point c3 = cubic_bezier[3] * t;

40742313779ee5e43be93a9191f1c86412cf183bKrzysztof Kosiński cubic_bbox(c0[0], c0[1], c1[0], c1[1], c2[0], c2[1], c3[0], c3[1], bbox);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen } else {

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen // should handle all not-so-easy curves:

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen Geom::Curve *ctemp = cit->transformed(t);

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen bbox.unionWith( ctemp->boundsExact());

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen delete ctemp;

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen }

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen }

40742313779ee5e43be93a9191f1c86412cf183bKrzysztof Kosiński }

981b809bc6ed10a21e89444d9447e5475801874fjohanengelen //return Geom::bounds_exact(pv * t);

40742313779ee5e43be93a9191f1c86412cf183bKrzysztof Kosiński return bbox;

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen}

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelenstatic void

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelengeom_line_wind_distance (Geom::Coord x0, Geom::Coord y0, Geom::Coord x1, Geom::Coord y1, Geom::Point const &pt, int *wind, Geom::Coord *best)

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen{

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen Geom::Coord Ax, Ay, Bx, By, Dx, Dy, s;

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen Geom::Coord dist2;

63267518b4ce196caab66ef8cbdcfc0921206b3djohanengelen

a4030d5ca449e7e384bc699cd249ee704faaeab0Chris Morgan /* Find distance */

Ax = x0;

Ay = y0;

Bx = x1;

By = y1;

Dx = x1 - x0;

Dy = y1 - y0;

const Geom::Coord Px = pt[X];

const Geom::Coord Py = pt[Y];

if (best) {

s = ((Px - Ax) * Dx + (Py - Ay) * Dy) / (Dx * Dx + Dy * Dy);

if (s <= 0.0) {

dist2 = (Px - Ax) * (Px - Ax) + (Py - Ay) * (Py - Ay);

} else if (s >= 1.0) {

dist2 = (Px - Bx) * (Px - Bx) + (Py - By) * (Py - By);

} else {

Geom::Coord Qx, Qy;

Qx = Ax + s * Dx;

Qy = Ay + s * Dy;

dist2 = (Px - Qx) * (Px - Qx) + (Py - Qy) * (Py - Qy);

}

if (dist2 < (*best * *best)) *best = sqrt (dist2);

}

if (wind) {

/* Find wind */

if ((Ax >= Px) && (Bx >= Px)) return;

if ((Ay >= Py) && (By >= Py)) return;

if ((Ay < Py) && (By < Py)) return;

if (Ay == By) return;

/* Ctach upper y bound */

if (Ay == Py) {

if (Ax < Px) *wind -= 1;

return;

} else if (By == Py) {

if (Bx < Px) *wind += 1;

return;

} else {

Geom::Coord Qx;

/* Have to calculate intersection */

Qx = Ax + Dx * (Py - Ay) / Dy;

if (Qx < Px) {

*wind += (Dy > 0.0) ? 1 : -1;

}

}

}

}

static void

geom_cubic_bbox_wind_distance (Geom::Coord x000, Geom::Coord y000,

Geom::Coord x001, Geom::Coord y001,

Geom::Coord x011, Geom::Coord y011,

Geom::Coord x111, Geom::Coord y111,

Geom::Point const &pt,

Geom::Rect *bbox, int *wind, Geom::Coord *best,

Geom::Coord tolerance)

{

Geom::Coord x0, y0, x1, y1, len2;

int needdist, needwind, needline;

const Geom::Coord Px = pt[X];

const Geom::Coord Py = pt[Y];

needdist = 0;

needwind = 0;

needline = 0;

if (bbox) cubic_bbox (x000, y000, x001, y001, x011, y011, x111, y111, *bbox);

x0 = std::min (x000, x001);

x0 = std::min (x0, x011);

x0 = std::min (x0, x111);

y0 = std::min (y000, y001);

y0 = std::min (y0, y011);

y0 = std::min (y0, y111);

x1 = std::max (x000, x001);

x1 = std::max (x1, x011);

x1 = std::max (x1, x111);

y1 = std::max (y000, y001);

y1 = std::max (y1, y011);

y1 = std::max (y1, y111);

if (best) {

/* Quickly adjust to endpoints */

len2 = (x000 - Px) * (x000 - Px) + (y000 - Py) * (y000 - Py);

if (len2 < (*best * *best)) *best = (Geom::Coord) sqrt (len2);

len2 = (x111 - Px) * (x111 - Px) + (y111 - Py) * (y111 - Py);

if (len2 < (*best * *best)) *best = (Geom::Coord) sqrt (len2);

if (((x0 - Px) < *best) && ((y0 - Py) < *best) && ((Px - x1) < *best) && ((Py - y1) < *best)) {

/* Point is inside sloppy bbox */

/* Now we have to decide, whether subdivide */

/* fixme: (Lauris) */

if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {

needdist = 1;

} else {

needline = 1;

}

}

}

if (!needdist && wind) {

if ((y1 >= Py) && (y0 < Py) && (x0 < Px)) {

/* Possible intersection at the left */

/* Now we have to decide, whether subdivide */

/* fixme: (Lauris) */

if (((y1 - y0) > 5.0) || ((x1 - x0) > 5.0)) {

needwind = 1;

} else {

needline = 1;

}

}

}

if (needdist || needwind) {

Geom::Coord x00t, x0tt, xttt, x1tt, x11t, x01t;

Geom::Coord y00t, y0tt, yttt, y1tt, y11t, y01t;

Geom::Coord s, t;

t = 0.5;

s = 1 - t;

x00t = s * x000 + t * x001;

x01t = s * x001 + t * x011;

x11t = s * x011 + t * x111;

x0tt = s * x00t + t * x01t;

x1tt = s * x01t + t * x11t;

xttt = s * x0tt + t * x1tt;

y00t = s * y000 + t * y001;

y01t = s * y001 + t * y011;

y11t = s * y011 + t * y111;

y0tt = s * y00t + t * y01t;

y1tt = s * y01t + t * y11t;

yttt = s * y0tt + t * y1tt;

geom_cubic_bbox_wind_distance (x000, y000, x00t, y00t, x0tt, y0tt, xttt, yttt, pt, NULL, wind, best, tolerance);

geom_cubic_bbox_wind_distance (xttt, yttt, x1tt, y1tt, x11t, y11t, x111, y111, pt, NULL, wind, best, tolerance);

} else if (1 || needline) {

geom_line_wind_distance (x000, y000, x111, y111, pt, wind, best);

}

}

static void

geom_curve_bbox_wind_distance(Geom::Curve const & c, Geom::Affine const &m,

Geom::Point const &pt,

Geom::Rect *bbox, int *wind, Geom::Coord *dist,

Geom::Coord tolerance, Geom::Rect const *viewbox,

Geom::Point &p0) // pass p0 through as it represents the last endpoint added (the finalPoint of last curve)

{

unsigned order = 0;

if (Geom::BezierCurve const* b = dynamic_cast<Geom::BezierCurve const*>(&c)) {

order = b->order();

}

if (order == 1) {

Geom::Point pe = c.finalPoint() * m;

if (bbox) {

bbox->expandTo(pe);

}

if (dist || wind) {

if (wind) { // we need to pick fill, so do what we're told

geom_line_wind_distance (p0[X], p0[Y], pe[X], pe[Y], pt, wind, dist);

} else { // only stroke is being picked; skip this segment if it's totally outside the viewbox

Geom::Rect swept(p0, pe);

if (!viewbox || swept.intersects(*viewbox))

geom_line_wind_distance (p0[X], p0[Y], pe[X], pe[Y], pt, wind, dist);

}

}

p0 = pe;

}

else if (order == 3) {

Geom::CubicBezier const& cubic_bezier = static_cast<Geom::CubicBezier const&>(c);

Geom::Point p1 = cubic_bezier[1] * m;

Geom::Point p2 = cubic_bezier[2] * m;

Geom::Point p3 = cubic_bezier[3] * m;

// get approximate bbox from handles (convex hull property of beziers):

Geom::Rect swept(p0, p3);

swept.expandTo(p1);

swept.expandTo(p2);

if (!viewbox || swept.intersects(*viewbox)) { // we see this segment, so do full processing

geom_cubic_bbox_wind_distance ( p0[X], p0[Y],

p1[X], p1[Y],

p2[X], p2[Y],

p3[X], p3[Y],

pt,

bbox, wind, dist, tolerance);

} else {

if (wind) { // if we need fill, we can just pretend it's a straight line

geom_line_wind_distance (p0[X], p0[Y], p3[X], p3[Y], pt, wind, dist);

} else { // otherwise, skip it completely

}

}

p0 = p3;

} else {

//this case handles sbasis as well as all other curve types

Geom::Path sbasis_path = Geom::cubicbezierpath_from_sbasis(c.toSBasis(), 0.1);

//recurse to convert the new path resulting from the sbasis to svgd

for (Geom::Path::iterator iter = sbasis_path.begin(); iter != sbasis_path.end(); ++iter) {

geom_curve_bbox_wind_distance(*iter, m, pt, bbox, wind, dist, tolerance, viewbox, p0);

}

}

}

pathv_matrix_point_bbox_wind_distance (Geom::PathVector const & pathv, Geom::Affine const &m, Geom::Point const &pt,

Geom::Rect *bbox, int *wind, Geom::Coord *dist,

Geom::Coord tolerance, Geom::Rect const *viewbox)

{

if (pathv.empty()) {

if (wind) *wind = 0;

if (dist) *dist = Geom::infinity();

return;

}

// remember last point of last curve

Geom::Point p0(0,0);

// remembering the start of subpath

Geom::Point p_start(0,0);

bool start_set = false;

for (Geom::PathVector::const_iterator it = pathv.begin(); it != pathv.end(); ++it) {

if (start_set) { // this is a new subpath

if (wind && (p0 != p_start)) // for correct fill picking, each subpath must be closed

geom_line_wind_distance (p0[X], p0[Y], p_start[X], p_start[Y], pt, wind, dist);

}

p0 = it->initialPoint() * m;

p_start = p0;

start_set = true;

if (bbox) {

bbox->expandTo(p0);

}

// loop including closing segment if path is closed

for (Geom::Path::const_iterator cit = it->begin(); cit != it->end_default(); ++cit) {

geom_curve_bbox_wind_distance(*cit, m, pt, bbox, wind, dist, tolerance, viewbox, p0);

}

}

if (start_set) {

if (wind && (p0 != p_start)) // for correct picking, each subpath must be closed

geom_line_wind_distance (p0[X], p0[Y], p_start[X], p_start[Y], pt, wind, dist);

}

}

Geom::PathVector

pathv_to_linear_and_cubic_beziers( Geom::PathVector const &pathv )

{

Geom::PathVector output;

for (Geom::PathVector::const_iterator pit = pathv.begin(); pit != pathv.end(); ++pit) {

output.push_back( Geom::Path() );

output.back().start( pit->initialPoint() );

output.back().close( pit->closed() );

for (Geom::Path::const_iterator cit = pit->begin(); cit != pit->end_open(); ++cit) {

if (is_straight_curve(*cit)) {

Geom::LineSegment l(cit->initialPoint(), cit->finalPoint());

output.back().append(l);

} else {

Geom::BezierCurve const *curve = dynamic_cast<Geom::BezierCurve const *>(&*cit);

if (curve && curve->order() == 3) {

Geom::CubicBezier b((*curve)[0], (*curve)[1], (*curve)[2], (*curve)[3]);

output.back().append(b);

} else {

// convert all other curve types to cubicbeziers

Geom::Path cubicbezier_path = Geom::cubicbezierpath_from_sbasis(cit->toSBasis(), 0.1);

output.back().append(cubicbezier_path);

}

}

}

}

return output;

}

Geom::PathVector

pathv_to_linear( Geom::PathVector const &pathv, double /*maxdisp*/)

{

Geom::PathVector output;

Geom::PathVector tmppath = pathv_to_linear_and_cubic_beziers(pathv);

// Now all path segments are either already lines, or they are beziers.

for (Geom::PathVector::const_iterator pit = tmppath.begin(); pit != tmppath.end(); ++pit) {

output.push_back( Geom::Path() );

output.back().start( pit->initialPoint() );

output.back().close( pit->closed() );

for (Geom::Path::const_iterator cit = pit->begin(); cit != pit->end_open(); ++cit) {

if (is_straight_curve(*cit)) {

Geom::LineSegment ls(cit->initialPoint(), cit->finalPoint());

output.back().append(ls);

}

else { /* all others must be Bezier curves */

Geom::BezierCurve const *curve = dynamic_cast<Geom::BezierCurve const *>(&*cit);

Geom::CubicBezier b((*curve)[0], (*curve)[1], (*curve)[2], (*curve)[3]);

std::vector<Geom::Point> bzrpoints = b.points();

Geom::Point A = bzrpoints[0];

Geom::Point B = bzrpoints[1];

Geom::Point C = bzrpoints[2];

Geom::Point D = bzrpoints[3];

std::vector<Geom::Point> pointlist;

pointlist.push_back(A);

recursive_bezier4(

A[X], A[Y],

B[X], B[Y],

C[X], C[Y],

D[X], D[Y],

pointlist,

0);

pointlist.push_back(D);

Geom::Point r1 = pointlist[0];

for (unsigned int i=1; i<pointlist.size();i++){

Geom::Point prev_r1 = r1;

r1 = pointlist[i];

Geom::LineSegment ls(prev_r1, r1);

output.back().append(ls);

}

pointlist.clear();

}

}

}

return output;

}

recursive_bezier4(const double x1, const double y1,

const double x2, const double y2,

const double x3, const double y3,

const double x4, const double y4,

std::vector<Geom::Point> &m_points,

int level)

{

// some of these should be parameters, but do it this way for now.

const double curve_collinearity_epsilon = 1e-30;

const double curve_angle_tolerance_epsilon = 0.01;

double m_cusp_limit = 0.0;

double m_angle_tolerance = 0.0;

double m_approximation_scale = 1.0;

double m_distance_tolerance_square = 0.5 / m_approximation_scale;

m_distance_tolerance_square *= m_distance_tolerance_square;

enum curve_recursion_limit_e { curve_recursion_limit = 32 };

#define calc_sq_distance(A,B,C,D) ((A-C)*(A-C) + (B-D)*(B-D))

if(level > curve_recursion_limit)

{

return;

}

// Calculate all the mid-points of the line segments

//----------------------

double x12 = (x1 + x2) / 2;

double y12 = (y1 + y2) / 2;

double x23 = (x2 + x3) / 2;

double y23 = (y2 + y3) / 2;

double x34 = (x3 + x4) / 2;

double y34 = (y3 + y4) / 2;

double x123 = (x12 + x23) / 2;

double y123 = (y12 + y23) / 2;

double x234 = (x23 + x34) / 2;

double y234 = (y23 + y34) / 2;

double x1234 = (x123 + x234) / 2;

double y1234 = (y123 + y234) / 2;

// Try to approximate the full cubic curve by a single straight line

//------------------

double dx = x4-x1;

double dy = y4-y1;

double d2 = fabs(((x2 - x4) * dy - (y2 - y4) * dx));

double d3 = fabs(((x3 - x4) * dy - (y3 - y4) * dx));

double da1, da2, k;

switch((int(d2 > curve_collinearity_epsilon) << 1) +

int(d3 > curve_collinearity_epsilon))

{

case 0:

// All collinear OR p1==p4

//----------------------

k = dx*dx + dy*dy;

if(k == 0)

{

d2 = calc_sq_distance(x1, y1, x2, y2);

d3 = calc_sq_distance(x4, y4, x3, y3);

}

else

{

k = 1 / k;

da1 = x2 - x1;

da2 = y2 - y1;

d2 = k * (da1*dx + da2*dy);

da1 = x3 - x1;

da2 = y3 - y1;

d3 = k * (da1*dx + da2*dy);

if(d2 > 0 && d2 < 1 && d3 > 0 && d3 < 1)

{

// Simple collinear case, 1---2---3---4

// We can leave just two endpoints

return;

}

if(d2 <= 0) d2 = calc_sq_distance(x2, y2, x1, y1);

else if(d2 >= 1) d2 = calc_sq_distance(x2, y2, x4, y4);

else d2 = calc_sq_distance(x2, y2, x1 + d2*dx, y1 + d2*dy);

if(d3 <= 0) d3 = calc_sq_distance(x3, y3, x1, y1);

else if(d3 >= 1) d3 = calc_sq_distance(x3, y3, x4, y4);

else d3 = calc_sq_distance(x3, y3, x1 + d3*dx, y1 + d3*dy);

}

if(d2 > d3)

{

if(d2 < m_distance_tolerance_square)

{

m_points.push_back(Geom::Point(x2, y2));

return;

}

}

else

{

if(d3 < m_distance_tolerance_square)

{

m_points.push_back(Geom::Point(x3, y3));

return;

}

}

break;

case 1:

// p1,p2,p4 are collinear, p3 is significant

//----------------------

if(d3 * d3 <= m_distance_tolerance_square * (dx*dx + dy*dy))

{

if(m_angle_tolerance < curve_angle_tolerance_epsilon)

{

m_points.push_back(Geom::Point(x23, y23));

return;

}

// Angle Condition

//----------------------

da1 = fabs(atan2(y4 - y3, x4 - x3) - atan2(y3 - y2, x3 - x2));

if(da1 >= M_PI) da1 = 2*M_PI - da1;

if(da1 < m_angle_tolerance)

{

m_points.push_back(Geom::Point(x2, y2));

m_points.push_back(Geom::Point(x3, y3));

return;

}

if(m_cusp_limit != 0.0)

{

if(da1 > m_cusp_limit)

{

m_points.push_back(Geom::Point(x3, y3));

return;

}

}

}

break;

case 2:

// p1,p3,p4 are collinear, p2 is significant

//----------------------

if(d2 * d2 <= m_distance_tolerance_square * (dx*dx + dy*dy))

{

if(m_angle_tolerance < curve_angle_tolerance_epsilon)

{

m_points.push_back(Geom::Point(x23, y23));

return;

}

// Angle Condition

//----------------------

da1 = fabs(atan2(y3 - y2, x3 - x2) - atan2(y2 - y1, x2 - x1));

if(da1 >= M_PI) da1 = 2*M_PI - da1;

if(da1 < m_angle_tolerance)

{

m_points.push_back(Geom::Point(x2, y2));

m_points.push_back(Geom::Point(x3, y3));

return;

}

if(m_cusp_limit != 0.0)

{

if(da1 > m_cusp_limit)

{

m_points.push_back(Geom::Point(x2, y2));

return;

}

}

}

break;

case 3:

// Regular case

//-----------------

if((d2 + d3)*(d2 + d3) <= m_distance_tolerance_square * (dx*dx + dy*dy))

{

// If the curvature doesn't exceed the distance_tolerance value

// we tend to finish subdivisions.

//----------------------

if(m_angle_tolerance < curve_angle_tolerance_epsilon)

{

m_points.push_back(Geom::Point(x23, y23));

return;

}

// Angle & Cusp Condition

//----------------------

k = atan2(y3 - y2, x3 - x2);

da1 = fabs(k - atan2(y2 - y1, x2 - x1));

da2 = fabs(atan2(y4 - y3, x4 - x3) - k);

if(da1 >= M_PI) da1 = 2*M_PI - da1;

if(da2 >= M_PI) da2 = 2*M_PI - da2;

if(da1 + da2 < m_angle_tolerance)

{

// Finally we can stop the recursion

//----------------------

m_points.push_back(Geom::Point(x23, y23));

return;

}

if(m_cusp_limit != 0.0)

{

if(da1 > m_cusp_limit)

{

m_points.push_back(Geom::Point(x2, y2));

return;

}

if(da2 > m_cusp_limit)

{

m_points.push_back(Geom::Point(x3, y3));

return;

}

}

}

break;

}

// Continue subdivision

//----------------------

recursive_bezier4(x1, y1, x12, y12, x123, y123, x1234, y1234, m_points, level + 1);

recursive_bezier4(x1234, y1234, x234, y234, x34, y34, x4, y4, m_points, level + 1);

}

void round_rectangle_outwards(Geom::Rect & rect) {

Geom::Interval ints[2];

for (int i=0; i < 2; i++) {

ints[i] = Geom::Interval(std::floor(rect[i][0]), std::ceil(rect[i][1]));

}

rect = Geom::Rect(ints[0], ints[1]);

}

namespace Geom {

bool transform_equalp(Geom::Affine const &m0, Geom::Affine const &m1, Geom::Coord const epsilon) {

return

Geom::are_near(m0[0], m1[0], epsilon) &&

Geom::are_near(m0[1], m1[1], epsilon) &&

Geom::are_near(m0[2], m1[2], epsilon) &&

Geom::are_near(m0[3], m1[3], epsilon);

}

bool translate_equalp(Geom::Affine const &m0, Geom::Affine const &m1, Geom::Coord const epsilon) {

return Geom::are_near(m0[4], m1[4], epsilon) && Geom::are_near(m0[5], m1[5], epsilon);

}

bool matrix_equalp(Geom::Affine const &m0, Geom::Affine const &m1, Geom::Coord const epsilon) {

return transform_equalp(m0, m1, epsilon) && translate_equalp(m0, m1, epsilon);

}

} //end namespace Geom