1 /*
2 * Copyright 2012, 2013 Thomas Schöps
3 * Copyright 2012-2019 Kai Pastor
4 *
5 * This file is part of OpenOrienteering.
6 *
7 * OpenOrienteering is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * OpenOrienteering is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with OpenOrienteering. If not, see <http://www.gnu.org/licenses/>.
19 */
20
21
22 #include "object.h"
23
24 #include <algorithm>
25 #include <cmath>
26 #include <cstddef>
27 #include <iterator>
28 #include <memory>
29 #include <stdexcept>
30 #include <type_traits>
31
32 #include <Qt>
33 #include <QtMath>
34 #include <QtNumeric>
35 #include <QFlags>
36 #include <QLatin1String>
37 #include <QPoint>
38 #include <QPointF>
39 #include <QScopedPointer>
40 #include <QStringRef>
41 #include <QTransform>
42 #include <QXmlStreamReader>
43 #include <QXmlStreamWriter>
44
45 #include <private/qbezier_p.h>
46
47 #include "settings.h"
48 #include "core/map.h"
49 #include "core/objects/text_object.h"
50 #include "core/renderables/renderable.h"
51 #include "core/symbols/line_symbol.h"
52 #include "core/symbols/point_symbol.h"
53 #include "core/symbols/symbol.h"
54 #include "core/symbols/text_symbol.h"
55 #include "core/virtual_coord_vector.h"
56 #include "fileformats/file_format.h"
57 #include "fileformats/file_import_export.h"
58 #include "util/util.h"
59 #include "util/xml_stream_util.h"
60
61 class QRectF;
62
63
64 // ### A namespace which collects various string constants of type QLatin1String. ###
65
66 namespace literal
67 {
68 static const QLatin1String object("object");
69 static const QLatin1String symbol("symbol");
70 static const QLatin1String type("type");
71 static const QLatin1String text("text");
72 static const QLatin1String coords("coords");
73 static const QLatin1String h_align("h_align");
74 static const QLatin1String v_align("v_align");
75 static const QLatin1String pattern("pattern");
76 static const QLatin1String rotation("rotation");
77 static const QLatin1String size("size");
78 static const QLatin1String tags("tags");
79 }
80
81
82
83 namespace OpenOrienteering {
84
85 // ### Object implementation ###
86
Object(Object::Type type,const Symbol * symbol)87 Object::Object(Object::Type type, const Symbol* symbol)
88 : type(type)
89 , symbol(symbol)
90 , output(*this)
91 {
92 // nothing
93 }
94
Object(Object::Type type,const Symbol * symbol,MapCoordVector coords,Map * map)95 Object::Object(Object::Type type, const Symbol* symbol, MapCoordVector coords, Map* map)
96 : type(type)
97 , symbol(symbol)
98 , coords(std::move(coords))
99 , map(map)
100 , output(*this)
101 {
102 // nothing
103 }
104
Object(const Object & proto)105 Object::Object(const Object& proto)
106 : type(proto.type)
107 , symbol(proto.symbol)
108 , coords(proto.coords)
109 , object_tags(proto.object_tags)
110 , rotation(proto.rotation)
111 , extent(proto.extent)
112 , output(*this)
113 {
114 // nothing
115 }
116
117 Object::~Object() = default;
118
copyFrom(const Object & other)119 void Object::copyFrom(const Object& other)
120 {
121 if (&other == this)
122 return;
123
124 if (type != other.type)
125 throw std::invalid_argument(Q_FUNC_INFO);
126
127 symbol = other.symbol;
128 coords = other.coords;
129 rotation = other.rotation;
130 // map unchanged!
131 object_tags = other.object_tags;
132 output_dirty = true;
133 extent = other.extent;
134 }
135
equals(const Object * other,bool compare_symbol) const136 bool Object::equals(const Object* other, bool compare_symbol) const
137 {
138 if (type != other->type)
139 return false;
140 if (compare_symbol)
141 {
142 if (bool(symbol) != bool(other->symbol))
143 return false;
144 if (symbol && !symbol->equals(other->symbol))
145 return false;
146 }
147
148 if (type == Text)
149 {
150 if (coords.front() != other->coords.front())
151 return false;
152 }
153 else
154 {
155 if (coords.size() != other->coords.size())
156 return false;
157 for (size_t i = 0, end = coords.size(); i < end; ++i)
158 {
159 if (coords[i] != other->coords[i])
160 return false;
161 }
162 }
163
164 if (object_tags != other->object_tags)
165 return false;
166
167 if (qAbs(getRotation() - other->getRotation()) >= 0.000001) // six decimal places in XML
168 return false;
169
170 if (type == Path)
171 {
172 auto const* path_this = static_cast<PathObject const*>(this);
173 auto const* path_other = static_cast<PathObject const*>(other);
174 if (path_this->getPatternOrigin() != path_other->getPatternOrigin())
175 return false;
176 }
177 else if (type == Text)
178 {
179 auto const* text_this = static_cast<TextObject const*>(this);
180 auto const* text_other = static_cast<TextObject const*>(other);
181
182 if (text_this->getBoxSize() != text_other->getBoxSize())
183 return false;
184 if (text_this->getText().compare(text_other->getText(), Qt::CaseSensitive) != 0)
185 return false;
186 if (text_this->getHorizontalAlignment() != text_other->getHorizontalAlignment())
187 return false;
188 if (text_this->getVerticalAlignment() != text_other->getVerticalAlignment())
189 return false;
190 }
191
192 return true;
193 }
194
195
196
validate() const197 bool Object::validate() const
198 {
199 return true;
200 }
201
202
203
asPoint()204 PointObject* Object::asPoint()
205 {
206 Q_ASSERT(type == Point);
207 return static_cast<PointObject*>(this);
208 }
209
asPoint() const210 const PointObject* Object::asPoint() const
211 {
212 Q_ASSERT(type == Point);
213 return static_cast<PointObject const*>(this);
214 }
215
asPath()216 PathObject* Object::asPath()
217 {
218 Q_ASSERT(type == Path);
219 return static_cast<PathObject*>(this);
220 }
221
asPath() const222 const PathObject* Object::asPath() const
223 {
224 Q_ASSERT(type == Path);
225 return static_cast<PathObject const*>(this);
226 }
227
asText()228 TextObject* Object::asText()
229 {
230 Q_ASSERT(type == Text);
231 return static_cast<TextObject*>(this);
232 }
233
asText() const234 const TextObject* Object::asText() const
235 {
236 Q_ASSERT(type == Text);
237 return static_cast<TextObject const*>(this);
238 }
239
240
241
save(QXmlStreamWriter & xml) const242 void Object::save(QXmlStreamWriter& xml) const
243 {
244 XmlElementWriter object_element(xml, literal::object);
245 object_element.writeAttribute(literal::type, type);
246 int symbol_index = -1;
247 if (map)
248 symbol_index = map->findSymbolIndex(symbol);
249 if (symbol_index != -1)
250 object_element.writeAttribute(literal::symbol, symbol_index);
251 if (symbol->isRotatable() && !qIsNull(rotation))
252 object_element.writeAttribute(literal::rotation, rotation);
253
254 if (type == Text)
255 {
256 auto const* text = static_cast<TextObject const*>(this);
257 object_element.writeAttribute(literal::h_align, text->getHorizontalAlignment());
258 object_element.writeAttribute(literal::v_align, text->getVerticalAlignment());
259 // For compatibility, we must keep the box size in the second coord ATM.
260 /// \todo From Mapper 1.0, save just the single anchor coordinate.
261 auto* object = const_cast<Object*>(this);
262 if (text->hasSingleAnchor())
263 {
264 object->coords.resize(1);
265 }
266 else
267 {
268 object->coords.resize(2);
269 object->coords.back() = text->getBoxSize();
270 }
271 }
272
273 if (!object_tags.empty())
274 {
275 XmlElementWriter tags_element(xml, literal::tags);
276 tags_element.write(object_tags);
277 }
278
279 {
280 // Scope of coords XML element
281 XmlElementWriter coords_element(xml, literal::coords);
282 coords_element.write(coords);
283 }
284
285 if (type == Path)
286 {
287 auto const* path = static_cast<PathObject const*>(this);
288 XmlElementWriter pattern_element(xml, literal::pattern);
289 pattern_element.writeAttribute(literal::rotation, path->getPatternRotation());
290 path->getPatternOrigin().save(xml);
291 }
292 else if (type == Text)
293 {
294 auto const* text = static_cast<TextObject const*>(this);
295 if (!text->hasSingleAnchor())
296 {
297 XmlElementWriter size_element(xml, literal::size);
298 auto size = text->getBoxSize();
299 size_element.writeAttribute(XmlStreamLiteral::width, size.nativeX());
300 size_element.writeAttribute(XmlStreamLiteral::height, size.nativeY());
301 }
302 xml.writeTextElement(literal::text, text->getText());
303 }
304 }
305
load(QXmlStreamReader & xml,Map * map,const SymbolDictionary & symbol_dict,const Symbol * symbol)306 Object* Object::load(QXmlStreamReader& xml, Map* map, const SymbolDictionary& symbol_dict, const Symbol* symbol)
307 {
308 Q_ASSERT(xml.name() == literal::object);
309
310 XmlElementReader object_element(xml);
311
312 Object::Type object_type = object_element.attribute<Object::Type>(literal::type);
313 Object* object = Object::getObjectForType(object_type);
314 if (!object)
315 throw FileFormatException(::OpenOrienteering::ImportExport::tr("Error while loading an object of type %1.").arg(object_type));
316
317 object->map = map;
318
319 if (symbol)
320 object->symbol = symbol;
321 else
322 {
323 QString symbol_id = object_element.attribute<QString>(literal::symbol);
324 bool conversion_ok;
325 const auto id_converted = symbol_id.toInt(&conversion_ok);
326 if (!symbol_id.isEmpty() && !conversion_ok)
327 throw FileFormatException(::OpenOrienteering::ImportExport::tr("Malformed symbol ID '%1' at line %2 column %3.")
328 .arg(symbol_id).arg(xml.lineNumber())
329 .arg(xml.columnNumber()));
330 object->symbol = symbol_dict[id_converted]; // FIXME: cannot work for forward references
331 // NOTE: object->symbol may be nullptr.
332 }
333
334 if (!object->symbol || !object->symbol->isTypeCompatibleTo(object))
335 {
336 // Throwing an exception will cause loading to fail.
337 // Rather than not loading the broken file at all,
338 // use the same symbols which are used for importing GPS tracks etc.
339 // FIXME: Implement a way to send a warning to the user.
340 switch (object_type)
341 {
342 case Point:
343 object->symbol = Map::getUndefinedPoint();
344 break;
345 case Path:
346 object->symbol = Map::getUndefinedLine();
347 break;
348 case Text:
349 object->symbol = Map::getUndefinedText();
350 break;
351 default:
352 throw FileFormatException(
353 ::OpenOrienteering::ImportExport::tr("Unable to find symbol for object at %1:%2.").
354 arg(xml.lineNumber()).arg(xml.columnNumber()) );
355 }
356 }
357
358 if (object->symbol->isRotatable())
359 {
360 object->rotation = object_element.attribute<qreal>(literal::rotation);
361 }
362
363 if (object_type == Text)
364 {
365 auto* text = static_cast<TextObject*>(object);
366 text->setHorizontalAlignment(object_element.attribute<TextObject::HorizontalAlignment>(literal::h_align));
367 text->setVerticalAlignment(object_element.attribute<TextObject::VerticalAlignment>(literal::v_align));
368 }
369
370 while (xml.readNextStartElement())
371 {
372 if (xml.name() == literal::coords)
373 {
374 XmlElementReader coords_element(xml);
375 try {
376 if (object_type == Text)
377 {
378 coords_element.readForText(object->coords);
379 if (object->coords.size() > 1)
380 static_cast<TextObject*>(object)->setBoxSize(object->coords[1]);
381 }
382 else
383 {
384 coords_element.read(object->coords);
385 }
386 }
387 catch (FileFormatException& e)
388 {
389 throw FileFormatException(::OpenOrienteering::ImportExport::tr("Error while loading an object of type %1 at %2:%3: %4").
390 arg(object_type).arg(xml.lineNumber()).arg(xml.columnNumber()).arg(e.message()));
391 }
392 }
393 else if (xml.name() == literal::pattern && object_type == Path)
394 {
395 XmlElementReader element(xml);
396
397 auto* path = static_cast<PathObject*>(object);
398 path->setPatternRotation(element.attribute<qreal>(literal::rotation));
399 while (xml.readNextStartElement())
400 {
401 if (xml.name() == XmlStreamLiteral::coord)
402 {
403 try
404 {
405 path->setPatternOrigin(MapCoord::load(xml));
406 }
407 catch (const std::range_error& e)
408 {
409 /// \todo Add a warning, but don't throw - throwing lets loading fail.
410 // throw FileFormatException(::OpenOrienteering::MapCoord::tr(e.what()));
411 qDebug("%s", e.what());
412 }
413 }
414 else
415 {
416 xml.skipCurrentElement(); // unknown
417 }
418 }
419 }
420 else if (xml.name() == literal::text && object_type == Text)
421 {
422 auto* text = static_cast<TextObject*>(object);
423 text->setText(xml.readElementText());
424 }
425 else if (xml.name() == literal::size && object_type == Text)
426 {
427 auto* text = static_cast<TextObject*>(object);
428 XmlElementReader size_element(xml);
429 auto const w = size_element.attribute<decltype(MapCoord().nativeX())>(XmlStreamLiteral::width);
430 auto const h = size_element.attribute<decltype(MapCoord().nativeY())>(XmlStreamLiteral::height);
431 text->setBoxSize(MapCoord::fromNative(w, h));
432 }
433 else if (xml.name() == literal::tags)
434 {
435 XmlElementReader(xml).read(object->object_tags);
436 }
437 else
438 xml.skipCurrentElement(); // unknown
439 }
440
441 if (object_type == Path)
442 {
443 auto* path = static_cast<PathObject*>(object);
444 path->recalculateParts();
445 }
446 object->output_dirty = true;
447
448 if (map &&
449 ( object->coords.empty()
450 || !object->coords.front().isRegular()
451 || !object->coords.back().isRegular() ) )
452 {
453 map->markAsIrregular(object);
454 }
455
456 return object;
457 }
458
459
setRotation(qreal new_rotation)460 void Object::setRotation(qreal new_rotation)
461 {
462 if (!qIsNaN(new_rotation))
463 {
464 rotation = new_rotation;
465 setOutputDirty();
466 }
467 }
468
469
forceUpdate() const470 void Object::forceUpdate() const
471 {
472 output_dirty = true;
473 update();
474 }
475
update() const476 bool Object::update() const
477 {
478 if (!output_dirty)
479 return false;
480
481 Symbol::RenderableOptions options = Symbol::RenderNormal;
482 if (map)
483 {
484 options = QFlag(map->renderableOptions());
485 if (extent.isValid())
486 map->setObjectAreaDirty(extent);
487 }
488
489 output.deleteRenderables();
490
491 extent = QRectF();
492
493 updateEvent();
494
495 createRenderables(output, options);
496
497 Q_ASSERT(extent.right() < 60000000); // assert if bogus values are returned
498 output_dirty = false;
499
500 if (map)
501 {
502 map->insertRenderablesOfObject(this);
503 if (extent.isValid())
504 map->setObjectAreaDirty(extent);
505 }
506
507 return true;
508 }
509
updateEvent() const510 void Object::updateEvent() const
511 {
512 // nothing here
513 }
514
createRenderables(ObjectRenderables & output,Symbol::RenderableOptions options) const515 void Object::createRenderables(ObjectRenderables& output, Symbol::RenderableOptions options) const
516 {
517 symbol->createRenderables(this, VirtualCoordVector(coords), output, options);
518 }
519
move(qint32 dx,qint32 dy)520 void Object::move(qint32 dx, qint32 dy)
521 {
522 for (MapCoord& coord : coords)
523 {
524 coord.setNativeX(dx + coord.nativeX());
525 coord.setNativeY(dy + coord.nativeY());
526 }
527
528 setOutputDirty();
529 }
530
move(const MapCoord & offset)531 void Object::move(const MapCoord& offset)
532 {
533 for (MapCoord& coord : coords)
534 {
535 coord += offset;
536 }
537
538 setOutputDirty();
539 }
540
scale(const MapCoordF & center,double factor)541 void Object::scale(const MapCoordF& center, double factor)
542 {
543 for (MapCoord& coord : coords)
544 {
545 coord.setX(center.x() + (coord.x() - center.x()) * factor);
546 coord.setY(center.y() + (coord.y() - center.y()) * factor);
547 }
548
549 setOutputDirty();
550 }
551
scale(double factor_x,double factor_y)552 void Object::scale(double factor_x, double factor_y)
553 {
554 for (MapCoord& coord : coords)
555 {
556 coord.setX(coord.x() * factor_x);
557 coord.setY(coord.y() * factor_y);
558 }
559
560 setOutputDirty();
561 }
562
rotateAround(const MapCoordF & center,qreal angle)563 void Object::rotateAround(const MapCoordF& center, qreal angle)
564 {
565 auto sin_angle = std::sin(angle);
566 auto cos_angle = std::cos(angle);
567
568 for (auto& coord : coords)
569 {
570 auto const center_to_coord = MapCoordF(coord.x() - center.x(), coord.y() - center.y());
571 coord.setX(center.x() + cos_angle * center_to_coord.x() + sin_angle * center_to_coord.y());
572 coord.setY(center.y() - sin_angle * center_to_coord.x() + cos_angle * center_to_coord.y());
573 }
574
575 if (symbol->isRotatable())
576 {
577 rotation += angle;
578 }
579 setOutputDirty();
580 }
581
rotate(qreal angle)582 void Object::rotate(qreal angle)
583 {
584 auto sin_angle = std::sin(angle);
585 auto cos_angle = std::cos(angle);
586
587 for (auto& coord : coords)
588 {
589 auto const old_coord = MapCoordF(coord);
590 coord.setX(+ cos_angle * old_coord.x() + sin_angle * old_coord.y());
591 coord.setY(- sin_angle * old_coord.x() + cos_angle * old_coord.y());
592 }
593
594 if (symbol->isRotatable())
595 {
596 rotation += angle;
597 }
598 setOutputDirty();
599 }
600
601
isPointOnObject(const MapCoordF & coord,qreal tolerance,bool treat_areas_as_paths,bool extended_selection) const602 int Object::isPointOnObject(const MapCoordF& coord, qreal tolerance, bool treat_areas_as_paths, bool extended_selection) const
603 {
604 Symbol::Type type = symbol->getType();
605 auto contained_types = symbol->getContainedTypes();
606
607 // Points
608 if (type == Symbol::Point)
609 {
610 if (extended_selection)
611 return extent.contains(coord) ? Symbol::Point : Symbol::NoSymbol;
612
613 return (coord.distanceSquaredTo(MapCoordF(coords[0])) <= tolerance) ? Symbol::Point : Symbol::NoSymbol;
614 }
615
616 // First check using extent
617 auto extent_extension = ((contained_types & Symbol::Line) || treat_areas_as_paths) ? tolerance : 0;
618 if (coord.x() < extent.left() - extent_extension) return Symbol::NoSymbol;
619 if (coord.y() < extent.top() - extent_extension) return Symbol::NoSymbol;
620 if (coord.x() > extent.right() + extent_extension) return Symbol::NoSymbol;
621 if (coord.y() > extent.bottom() + extent_extension) return Symbol::NoSymbol;
622
623 if (type == Symbol::Text)
624 {
625 // Texts
626 auto const* text_object = static_cast<TextObject const*>(this);
627 return (text_object->calcTextPositionAt(coord, true) != -1) ? Symbol::Text : Symbol::NoSymbol;
628 }
629
630 // Path objects
631 auto const* path = static_cast<PathObject const*>(this);
632 return path->isPointOnPath(coord, tolerance, treat_areas_as_paths, true);
633 }
634
takeRenderables()635 void Object::takeRenderables()
636 {
637 output.takeRenderables();
638 }
639
clearRenderables()640 void Object::clearRenderables()
641 {
642 output.deleteRenderables();
643 extent = QRectF();
644 }
645
setSymbol(const Symbol * new_symbol,bool no_checks)646 bool Object::setSymbol(const Symbol* new_symbol, bool no_checks)
647 {
648 if (!no_checks && new_symbol)
649 {
650 if (!new_symbol->isTypeCompatibleTo(this))
651 return false;
652 }
653
654 symbol = new_symbol;
655 setOutputDirty();
656 return true;
657 }
658
getObjectForType(Object::Type type,const Symbol * symbol)659 Object* Object::getObjectForType(Object::Type type, const Symbol* symbol)
660 {
661 switch(type)
662 {
663 case Point:
664 return new PointObject(symbol);
665 case Path:
666 return new PathObject(symbol);
667 case Text:
668 return new TextObject(symbol);
669 }
670 return nullptr;
671 }
672
setTags(const Object::Tags & tags)673 void Object::setTags(const Object::Tags& tags)
674 {
675 if (object_tags != tags)
676 {
677 object_tags = tags;
678 if (map)
679 {
680 map->setObjectsDirty();
681 if (map->isObjectSelected(this))
682 map->emitSelectionEdited();
683 }
684 }
685 }
686
setTag(const QString & key,const QString & value)687 void Object::setTag(const QString& key, const QString& value)
688 {
689 if (!object_tags.contains(key) || object_tags.value(key) != value)
690 {
691 object_tags.insert(key, value);
692 if (map)
693 {
694 map->setObjectsDirty();
695 if (map->isObjectSelected(this))
696 map->emitSelectionEdited();
697 }
698 }
699 }
700
removeTag(const QString & key)701 void Object::removeTag(const QString& key)
702 {
703 if (object_tags.contains(key))
704 {
705 object_tags.remove(key);
706 if (map)
707 map->setObjectsDirty();
708 }
709 }
710
711
includeControlPointsRect(QRectF & rect) const712 void Object::includeControlPointsRect(QRectF& rect) const
713 {
714 if (type == Object::Path)
715 {
716 const PathObject* path = asPath();
717 for (auto const& coord : path->coords)
718 rectInclude(rect, QPointF(coord));
719 }
720 else if (type == Object::Text)
721 {
722 const TextObject* text = asText();
723 std::vector<QPointF> text_handles(text->controlPoints());
724 for (auto& text_handle : text_handles)
725 rectInclude(rect, text_handle);
726 }
727 }
728
729
730
731 // ### PathPart ###
732
733 /*
734 * Some headers may not want to include object.h but rather rely on
735 * PathPartVector::size_type being a std::size_t.
736 */
737 static_assert(std::is_same<PathPartVector::size_type, std::size_t>::value,
738 "PathCoordVector::size_type is not std::size_t");
739
PathPart(PathObject & path,const VirtualPath & virtual_path)740 PathPart::PathPart(PathObject& path, const VirtualPath& virtual_path)
741 : PathPart(path, virtual_path.first_index, virtual_path.last_index)
742 {
743 if (!virtual_path.path_coords.empty())
744 {
745 path_coords.reserve(virtual_path.path_coords.size());
746 path_coords.insert(end(path_coords), begin(virtual_path.path_coords), end(virtual_path.path_coords));
747 }
748 }
749
setClosed(bool closed,bool may_use_existing_close_point)750 void PathPart::setClosed(bool closed, bool may_use_existing_close_point)
751 {
752 Q_ASSERT(!empty());
753
754 if (!isClosed() && closed)
755 {
756 if (size() == 1 || !may_use_existing_close_point ||
757 !path->coords[first_index].isPositionEqualTo(path->coords[last_index]))
758 {
759 path->coords[last_index].setHolePoint(false);
760 path->coords[last_index].setClosePoint(false);
761 path->coords.insert(path->coords.begin() + last_index + 1, path->coords[first_index]);
762 path->partSizeChanged(path->findPartForIndex(first_index), 1);
763 }
764 path->setClosingPoint(last_index, path->coords[first_index]);
765 path->setOutputDirty();
766 }
767 else if (isClosed() && !closed)
768 {
769 path->coords[last_index].setClosePoint(false);
770 path->setOutputDirty();
771 }
772 }
773
connectEnds()774 void PathPart::connectEnds()
775 {
776 if (!isClosed())
777 {
778 path->coords[first_index] += path->coords[last_index];
779 path->coords[first_index] /= 2;
780 path->setClosingPoint(last_index, path->coords[first_index]);
781 path->setOutputDirty();
782 }
783 }
784
reverse()785 void PathPart::reverse()
786 {
787 MapCoordVector& coords = path->coords;
788
789 bool set_last_hole_point = false;
790 auto half = (first_index + last_index + 1) / 2;
791 for (auto c = first_index; c <= last_index; ++c)
792 {
793 MapCoord coord = coords[c];
794 if (c < half)
795 {
796 auto mirror_index = last_index - (c - first_index);
797 coords[c] = coords[mirror_index];
798 coords[mirror_index] = coord;
799 }
800
801 if (!(c == first_index && isClosed()) && coords[c].isCurveStart())
802 {
803 Q_ASSERT((c - first_index) >= 3);
804 coords[c - 3].setCurveStart(true);
805 if (!(c == last_index && isClosed()))
806 coords[c].setCurveStart(false);
807 }
808 else if (coords[c].isHolePoint())
809 {
810 coords[c].setHolePoint(false);
811 if (c >= first_index + 1)
812 coords[c - 1].setHolePoint(true);
813 else
814 set_last_hole_point = true;
815 }
816 }
817
818 if (coords[first_index].isClosePoint())
819 {
820 coords[first_index].setClosePoint(false);
821 coords[last_index].setClosePoint(true);
822 }
823
824 if (set_last_hole_point)
825 coords[last_index].setHolePoint(true);
826
827 path->setOutputDirty();
828 }
829
830 // static
calculatePathParts(const VirtualCoordVector & coords)831 PathPartVector PathPart::calculatePathParts(const VirtualCoordVector& coords)
832 {
833 PathPartVector parts;
834 PathCoordVector path_coords(coords);
835 auto last = coords.size();
836 auto part_start = VirtualCoordVector::size_type { 0 };
837 while (part_start != last)
838 {
839 auto part_end = path_coords.update(part_start);
840 parts.emplace_back(coords, part_start, part_end);
841 parts.back().path_coords.swap(path_coords);
842 part_start = part_end + 1;
843 }
844 return parts;
845 }
846
847
848
849 //### PathPartVector ###
850
851 // Not inline because it is to be used by reference.
compareEndIndex(const PathPart & part,VirtualPath::size_type index)852 bool PathPartVector::compareEndIndex(const PathPart& part, VirtualPath::size_type index)
853 {
854 return part.last_index+1 <= index;
855 }
856
857
858
859 // ### PathObject ###
860
PathObject(const Symbol * symbol)861 PathObject::PathObject(const Symbol* symbol)
862 : Object(Object::Path, symbol)
863 {
864 Q_ASSERT(!symbol || (symbol->getType() == Symbol::Line || symbol->getType() == Symbol::Area || symbol->getType() == Symbol::Combined));
865 }
866
PathObject(const Symbol * symbol,MapCoordVector coords,Map * map)867 PathObject::PathObject(const Symbol* symbol, MapCoordVector coords, Map* map)
868 : Object(Object::Path, symbol, std::move(coords), map)
869 {
870 Q_ASSERT(!symbol || (symbol->getType() == Symbol::Line || symbol->getType() == Symbol::Area || symbol->getType() == Symbol::Combined));
871 recalculateParts();
872 }
873
PathObject(const Symbol * symbol,const PathObject & proto,MapCoordVector::size_type piece)874 PathObject::PathObject(const Symbol* symbol, const PathObject& proto, MapCoordVector::size_type piece)
875 : Object { Object::Path, symbol }
876 {
877 auto begin = proto.coords.begin() + piece;
878 auto part = proto.findPartForIndex(piece);
879 if (piece == part->last_index)
880 {
881 if (part->isClosed())
882 begin = proto.coords.begin() + part->first_index;
883 else
884 begin = proto.coords.begin() + part->prevCoordIndex(piece);
885 }
886
887 auto end = begin + (begin->isCurveStart() ? 4 : 2);
888
889 std::for_each(begin, end, [this](MapCoord c)
__anonc09b97f00102(MapCoord c) 890 {
891 c.setHolePoint(false);
892 c.setClosePoint(false);
893 addCoordinate(c);
894 });
895
896 coords.back().setHolePoint(true);
897 }
898
PathObject(const PathObject & proto)899 PathObject::PathObject(const PathObject& proto)
900 : Object(proto)
901 , pattern_origin(proto.pattern_origin)
902 {
903 path_parts.reserve(proto.path_parts.size());
904 for (const PathPart& part : proto.path_parts)
905 {
906 path_parts.emplace_back(*this, part);
907 }
908 }
909
PathObject(const PathPart & proto_part)910 PathObject::PathObject(const PathPart &proto_part)
911 : Object(*proto_part.path)
912 , pattern_origin(proto_part.path->pattern_origin)
913 {
914 auto begin = proto_part.path->coords.begin();
915 coords.reserve(proto_part.size());
916 coords.assign(begin + proto_part.first_index, begin + (proto_part.last_index+1));
917 path_parts.emplace_back(*this, proto_part);
918 path_parts.front().last_index -= path_parts.front().first_index;
919 path_parts.front().first_index = 0;
920 }
921
duplicate() const922 PathObject* PathObject::duplicate() const
923 {
924 return new PathObject(*this);
925 }
926
copyFrom(const Object & other)927 void PathObject::copyFrom(const Object& other)
928 {
929 if (&other == this)
930 return;
931
932 Object::copyFrom(other);
933
934 const PathObject& other_path = *other.asPath();
935 pattern_origin = other_path.getPatternOrigin();
936
937 path_parts.clear();
938 path_parts.reserve(other_path.path_parts.size());
939 for (const PathPart& part : other_path.path_parts)
940 {
941 path_parts.emplace_back(*this, part);
942 }
943 }
944
945
946
validate() const947 bool PathObject::validate() const
948 {
949 return std::all_of(begin(path_parts), end(path_parts), [](auto &part) {
950 return !part.isClosed()
951 || part.coords[part.first_index] == part.coords[part.last_index];
952 });
953 }
954
955
956
normalize()957 void PathObject::normalize()
958 {
959 for (MapCoordVector::size_type i = 0; i < coords.size(); ++i)
960 {
961 if (coords[i].isCurveStart())
962 {
963 if (i+3 >= getCoordinateCount())
964 {
965 coords[i].setCurveStart(false);
966 continue;
967 }
968
969 if (coords[i + 1].isClosePoint() || coords[i + 1].isHolePoint() ||
970 coords[i + 2].isClosePoint() || coords[i + 2].isHolePoint())
971 {
972 coords[i].setCurveStart(false);
973 continue;
974 }
975
976 coords[i + 1].setCurveStart(false);
977 coords[i + 1].setDashPoint(false);
978 coords[i + 2].setCurveStart(false);
979 coords[i + 2].setDashPoint(false);
980 i += 2;
981 }
982
983 if (i > 0 && coords[i].isHolePoint())
984 {
985 if (coords[i-1].isHolePoint())
986 deleteCoordinate(i, false);
987 }
988 }
989 }
990
intersectsBox(const QRectF & box) const991 bool PathObject::intersectsBox(const QRectF& box) const
992 {
993 // Check path parts for an intersection with box
994 if (std::any_of(begin(path_parts), end(path_parts), [&box](const PathPart& part) { return part.intersectsBox(box); }))
995 {
996 return true;
997 }
998
999 // If this is an area, additionally check if the area contains the box
1000 if (getSymbol()->getContainedTypes() & Symbol::Area)
1001 {
1002 return isPointOnObject(MapCoordF(box.center()), 0, false, false);
1003 }
1004
1005 return false;
1006 }
1007
findPartForIndex(MapCoordVector::size_type coords_index) const1008 PathPartVector::const_iterator PathObject::findPartForIndex(MapCoordVector::size_type coords_index) const
1009 {
1010 return std::lower_bound(begin(path_parts), end(path_parts), coords_index, PathPartVector::compareEndIndex);
1011 }
1012
findPartForIndex(MapCoordVector::size_type coords_index)1013 PathPartVector::iterator PathObject::findPartForIndex(MapCoordVector::size_type coords_index)
1014 {
1015 setOutputDirty();
1016 return std::lower_bound(begin(path_parts), end(path_parts), coords_index, PathPartVector::compareEndIndex);
1017 }
1018
findPartIndexForIndex(MapCoordVector::size_type coords_index) const1019 PathPartVector::size_type PathObject::findPartIndexForIndex(MapCoordVector::size_type coords_index) const
1020 {
1021 for (PathPartVector::size_type i = 0; i < path_parts.size(); ++i)
1022 {
1023 if (path_parts[i].first_index <= coords_index && path_parts[i].last_index >= coords_index)
1024 return i;
1025 }
1026 Q_ASSERT(false);
1027 return 0;
1028 }
1029
findPathCoordForIndex(MapCoordVector::size_type index) const1030 PathCoord PathObject::findPathCoordForIndex(MapCoordVector::size_type index) const
1031 {
1032 auto part = findPartForIndex(index);
1033 if (part != end(path_parts))
1034 {
1035 return *(std::lower_bound(begin(part->path_coords), end(part->path_coords), index, PathCoord::indexLessThanValue));
1036 }
1037 return path_parts.front().path_coords.front();
1038 }
1039
isCurveHandle(MapCoordVector::size_type index) const1040 bool PathObject::isCurveHandle(MapCoordVector::size_type index) const
1041 {
1042 return ( index < coords.size() &&
1043 !coords[index].isCurveStart() &&
1044 (
1045 ( index > 0 && coords[index-1].isCurveStart() ) ||
1046 ( index > 1 && coords[index-2].isCurveStart() )
1047 )
1048 );
1049 }
1050
deletePart(PathPartVector::size_type part_index)1051 void PathObject::deletePart(PathPartVector::size_type part_index)
1052 {
1053 setOutputDirty();
1054
1055 auto part = begin(path_parts) + part_index;
1056 auto first_index = part->first_index;
1057 auto end_index = part->last_index + 1;
1058 auto first_coord = begin(coords);
1059 coords.erase(first_coord + first_index, first_coord + end_index);
1060 partSizeChanged(part, first_index - end_index);
1061 path_parts.erase(part);
1062 }
1063
partSizeChanged(PathPartVector::iterator part,MapCoordVector::difference_type change)1064 void PathObject::partSizeChanged(PathPartVector::iterator part, MapCoordVector::difference_type change)
1065 {
1066 Q_ASSERT(isOutputDirty());
1067
1068 part->last_index += change;
1069 auto last = end(path_parts);
1070 for (++part; part != last; ++part)
1071 {
1072 part->first_index += change;
1073 part->last_index += change;
1074 }
1075 }
1076
1077
transform(const QTransform & t)1078 void PathObject::transform(const QTransform& t)
1079 {
1080 if (t.isIdentity())
1081 return;
1082
1083 for (auto& coord : coords)
1084 {
1085 const auto p = t.map(MapCoordF{coord});
1086 coord.setX(p.x());
1087 coord.setY(p.y());
1088 }
1089 pattern_origin = MapCoord{t.map(MapCoordF{getPatternOrigin()})};
1090 setOutputDirty();
1091 }
1092
1093
1094
setPatternRotation(qreal rotation)1095 void PathObject::setPatternRotation(qreal rotation)
1096 {
1097 setRotation(rotation);
1098 }
1099
setPatternOrigin(const MapCoord & origin)1100 void PathObject::setPatternOrigin(const MapCoord& origin)
1101 {
1102 pattern_origin = origin;
1103 setOutputDirty();
1104 }
1105
calcClosestPointOnPath(MapCoordF coord,float & out_distance_sq,PathCoord & out_path_coord,MapCoordVector::size_type start_index,MapCoordVector::size_type end_index) const1106 void PathObject::calcClosestPointOnPath(
1107 MapCoordF coord,
1108 float& out_distance_sq,
1109 PathCoord& out_path_coord,
1110 MapCoordVector::size_type start_index,
1111 MapCoordVector::size_type end_index) const
1112 {
1113 update();
1114
1115 auto bound = std::numeric_limits<float>::max();
1116 for (const auto& part : path_parts)
1117 {
1118 if (part.first_index <= end_index && part.last_index >= start_index) /// \todo Legacy compatibility, review/remove
1119 {
1120 auto path_coord = part.findClosestPointTo(coord, out_distance_sq, bound, start_index, end_index);
1121 if (out_distance_sq < bound)
1122 {
1123 bound = out_distance_sq;
1124 out_path_coord = path_coord;
1125 }
1126 }
1127 }
1128 }
1129
calcClosestCoordinate(MapCoordF coord,float & out_distance_sq,MapCoordVector::size_type & out_index) const1130 void PathObject::calcClosestCoordinate(MapCoordF coord, float& out_distance_sq, MapCoordVector::size_type& out_index) const
1131 {
1132 update();
1133
1134 auto coords_size = coords.size();
1135 if (coords_size == 0)
1136 {
1137 out_distance_sq = -1;
1138 out_index = -1;
1139 return;
1140 }
1141
1142 // NOTE: do not try to optimize this by starting with index 1, it will overlook curve starts this way
1143 out_distance_sq = 999999;
1144 out_index = 0;
1145 for (MapCoordVector::size_type i = 0; i < coords_size; ++i)
1146 {
1147 double length_sq = (coord - MapCoordF(coords[i])).lengthSquared();
1148 if (length_sq < out_distance_sq)
1149 {
1150 out_distance_sq = length_sq;
1151 out_index = i;
1152 }
1153
1154 if (coords[i].isCurveStart())
1155 i += 2;
1156 }
1157 }
1158
subdivide(const PathCoord & path_coord)1159 MapCoordVector::size_type PathObject::subdivide(const PathCoord& path_coord)
1160 {
1161 return subdivide(path_coord.index, path_coord.param);
1162 }
1163
subdivide(MapCoordVector::size_type index,float param)1164 MapCoordVector::size_type PathObject::subdivide(MapCoordVector::size_type index, float param)
1165 {
1166 Q_ASSERT(index < coords.size());
1167
1168 if (coords[index].isCurveStart())
1169 {
1170 MapCoordF o0, o1, o2, o3, o4;
1171 PathCoord::splitBezierCurve(MapCoordF(coords[index]), MapCoordF(coords[index+1]),
1172 MapCoordF(coords[index+2]), MapCoordF(coords[index+3]),
1173 param, o0, o1, o2, o3, o4);
1174 coords[index + 1] = MapCoord(o0);
1175 coords[index + 2] = MapCoord(o4);
1176 addCoordinate(index + 2, MapCoord(o3));
1177 MapCoord middle_coord = MapCoord(o2);
1178 middle_coord.setCurveStart(true);
1179 addCoordinate(index + 2, middle_coord);
1180 addCoordinate(index + 2, MapCoord(o1));
1181 Q_ASSERT(isOutputDirty());
1182 return index + 3;
1183 }
1184
1185 addCoordinate(index + 1, MapCoord(MapCoordF(coords[index]) + (MapCoordF(coords[index+1]) - MapCoordF(coords[index])) * param));
1186 Q_ASSERT(isOutputDirty());
1187 return index + 1;
1188 }
1189
canBeConnected(const PathObject * other,double connect_threshold_sq) const1190 bool PathObject::canBeConnected(const PathObject* other, double connect_threshold_sq) const
1191 {
1192 for (const auto& part : path_parts)
1193 {
1194 if (part.isClosed())
1195 continue;
1196
1197 for (const auto& other_part : other->path_parts)
1198 {
1199 if (other_part.isClosed())
1200 continue;
1201
1202 if (coords[part.first_index].distanceSquaredTo(other->coords[other_part.first_index]) <= connect_threshold_sq
1203 || coords[part.first_index].distanceSquaredTo(other->coords[other_part.last_index]) <= connect_threshold_sq
1204 || coords[part.last_index].distanceSquaredTo(other->coords[other_part.first_index]) <= connect_threshold_sq
1205 || coords[part.last_index].distanceSquaredTo(other->coords[other_part.last_index]) <= connect_threshold_sq)
1206 {
1207 return true;
1208 }
1209 }
1210 }
1211
1212 return false;
1213 }
1214
connectIfClose(PathObject * other,double connect_threshold_sq)1215 bool PathObject::connectIfClose(PathObject* other, double connect_threshold_sq)
1216 {
1217 bool did_connect_path = false;
1218
1219 auto num_parts = parts().size();
1220 auto num_other_parts = other->parts().size();
1221 std::vector<bool> other_parts; // Which parts have not been connected to this part yet?
1222 other_parts.assign(num_other_parts, true);
1223 for (std::size_t i = 0; i < num_parts; ++i)
1224 {
1225 if (path_parts[i].isClosed())
1226 continue;
1227
1228 for (std::size_t k = 0; k < num_other_parts; ++k)
1229 {
1230 if (!other_parts[k] || other->path_parts[k].isClosed())
1231 continue;
1232
1233 if (coords[path_parts[i].first_index].distanceSquaredTo(other->coords[other->path_parts[k].first_index]) <= connect_threshold_sq)
1234 {
1235 other->path_parts[k].reverse();
1236 connectPathParts(i, other, k, true);
1237 }
1238 else if (coords[path_parts[i].first_index].distanceSquaredTo(other->coords[other->path_parts[k].last_index]) <= connect_threshold_sq)
1239 connectPathParts(i, other, k, true);
1240 else if (coords[path_parts[i].last_index].distanceSquaredTo(other->coords[other->path_parts[k].first_index]) <= connect_threshold_sq)
1241 connectPathParts(i, other, k, false);
1242 else if (coords[path_parts[i].last_index].distanceSquaredTo(other->coords[other->path_parts[k].last_index]) <= connect_threshold_sq)
1243 {
1244 other->path_parts[k].reverse();
1245 connectPathParts(i, other, k, false);
1246 }
1247 else
1248 continue;
1249
1250 if (coords[path_parts[i].first_index].distanceSquaredTo(coords[path_parts[i].last_index]) <= connect_threshold_sq)
1251 path_parts[i].connectEnds();
1252
1253 did_connect_path = true;
1254
1255 other_parts[k] = false;
1256 }
1257 }
1258
1259 if (did_connect_path)
1260 {
1261 // Copy over all remaining parts of the other object
1262 getCoordinateRef(getCoordinateCount() - 1).setHolePoint(true);
1263 for (std::size_t i = 0; i < num_other_parts; ++i)
1264 {
1265 if (other_parts[i])
1266 appendPathPart(other->parts()[i]);
1267 }
1268
1269 Q_ASSERT(isOutputDirty());
1270 }
1271
1272 return did_connect_path;
1273 }
1274
connectPathParts(PathPartVector::size_type part_index,const PathObject * other,PathPartVector::size_type other_part_index,bool prepend,bool merge_ends)1275 void PathObject::connectPathParts(PathPartVector::size_type part_index, const PathObject* other, PathPartVector::size_type other_part_index, bool prepend, bool merge_ends)
1276 {
1277 Q_ASSERT(part_index < path_parts.size());
1278 PathPart& part = path_parts[part_index];
1279 PathPart& other_part = other->path_parts[other_part_index];
1280 Q_ASSERT(!part.isClosed() && !other_part.isClosed());
1281
1282 auto const other_part_size = other_part.size();
1283 auto const appended_part_size = other_part_size - (merge_ends ? 1 : 0);
1284 coords.resize(coords.size() + appended_part_size);
1285
1286 if (prepend)
1287 {
1288 for (auto i = coords.size() - 1; i >= part.first_index + appended_part_size; --i)
1289 coords[i] = coords[i - appended_part_size];
1290
1291 MapCoord& join_coord = coords[part.first_index + appended_part_size];
1292 if (merge_ends)
1293 {
1294 join_coord.setNativeX((join_coord.nativeX() + other->coords[other_part.last_index].nativeX()) / 2);
1295 join_coord.setNativeY((join_coord.nativeY() + other->coords[other_part.last_index].nativeY()) / 2);
1296 }
1297 join_coord.setHolePoint(false);
1298 join_coord.setClosePoint(false);
1299
1300 for (auto i = part.first_index; i < part.first_index + appended_part_size; ++i)
1301 coords[i] = other->coords[i - part.first_index + other_part.first_index];
1302
1303 if (!merge_ends)
1304 {
1305 MapCoord& second_join_coord = coords[part.first_index + appended_part_size - 1];
1306 second_join_coord.setHolePoint(false);
1307 second_join_coord.setClosePoint(false);
1308 }
1309 }
1310 else
1311 {
1312 auto end_index = part.last_index;
1313
1314 if (merge_ends)
1315 {
1316 MapCoord coord = other->coords[other_part.first_index]; // take flags from first coord of path to append
1317 coord.setNativeX((coords[end_index].nativeX() + coord.nativeX()) / 2);
1318 coord.setNativeY((coords[end_index].nativeY() + coord.nativeY()) / 2);
1319 coords[end_index] = coord;
1320 }
1321 else
1322 {
1323 coords[end_index].setHolePoint(false);
1324 coords[end_index].setClosePoint(false);
1325 }
1326
1327 for (auto i = coords.size() - 1; i > part.last_index + appended_part_size; --i)
1328 coords[i] = coords[i - appended_part_size];
1329 for (auto i = part.last_index+1; i < part.last_index + 1 + appended_part_size; ++i)
1330 coords[i] = other->coords[i - end_index + other_part.first_index - (merge_ends ? 0 : 1)];
1331 }
1332
1333 setOutputDirty();
1334 partSizeChanged(begin(path_parts)+part_index, appended_part_size);
1335 Q_ASSERT(!part.isClosed());
1336 }
1337
removeFromLine(PathPartVector::size_type part_index,PathCoord::length_type removal_begin,PathCoord::length_type removal_end) const1338 std::vector<PathObject*> PathObject::removeFromLine(
1339 PathPartVector::size_type part_index,
1340 PathCoord::length_type removal_begin,
1341 PathCoord::length_type removal_end) const
1342 {
1343 Q_ASSERT(path_parts.size() == 1); // TODO
1344 Q_ASSERT((symbol->getContainedTypes() & ~Symbol::Combined) == Symbol::Line);
1345
1346 const PathPart& part = path_parts[part_index];
1347 Q_ASSERT(part.path_coords.front().clen <= removal_begin);
1348 Q_ASSERT(part.path_coords.front().clen <= removal_end);
1349 Q_ASSERT(part.path_coords.back().clen >= removal_begin);
1350 Q_ASSERT(part.path_coords.back().clen >= removal_end);
1351
1352 std::vector<PathObject*> objects;
1353
1354 if (part.path_coords.front().clen >= removal_begin
1355 && part.path_coords.back().clen <= removal_end)
1356 {
1357 // nothing left, or invalid
1358 }
1359 else if (removal_end < removal_begin || part.isClosed())
1360 {
1361 PathObject* obj = duplicate()->asPath();
1362 obj->changePathBounds(part_index, removal_end, removal_begin);
1363 objects.push_back(obj);
1364 }
1365 else
1366 {
1367 if (removal_begin > part.path_coords.front().clen)
1368 {
1369 PathObject* obj = duplicate()->asPath();
1370 obj->changePathBounds(part_index, part.path_coords.front().clen, removal_begin);
1371 objects.push_back(obj);
1372 }
1373 if (removal_end < part.path_coords.back().clen)
1374 {
1375 PathObject* obj = duplicate()->asPath();
1376 obj->changePathBounds(part_index, removal_end, part.path_coords.back().clen);
1377 objects.push_back(obj);
1378 }
1379 }
1380
1381 return objects;
1382 }
1383
splitLineAt(const PathCoord & split_pos) const1384 std::vector<PathObject*> PathObject::splitLineAt(const PathCoord& split_pos) const
1385 {
1386 Q_ASSERT(path_parts.size() == 1);
1387 Q_ASSERT((symbol->getContainedTypes() & ~Symbol::Combined) == Symbol::Line);
1388
1389 std::vector<PathObject*> objects;
1390 if (path_parts.size() == 1)
1391 {
1392 const auto part_index = PathPartVector::size_type { 0 };
1393 const auto& part = path_parts[part_index];
1394 if (part.isClosed())
1395 {
1396 PathObject* obj = duplicate()->asPath();
1397 if ( split_pos.clen == part.path_coords.front().clen ||
1398 split_pos.clen == part.path_coords.back().clen )
1399 {
1400 (obj->path_parts[part_index]).setClosed(false);
1401 }
1402 else
1403 {
1404 obj->changePathBounds(part_index, split_pos.clen, split_pos.clen);
1405 }
1406 objects.push_back(obj);
1407 }
1408 else if ( split_pos.clen != part.path_coords.front().clen &&
1409 split_pos.clen != part.path_coords.back().clen)
1410 {
1411 PathObject* obj1 = duplicate()->asPath();
1412 obj1->changePathBounds(part_index, part.path_coords.front().clen, split_pos.clen);
1413 objects.push_back(obj1);
1414
1415 PathObject* obj2 = duplicate()->asPath();
1416 obj2->changePathBounds(part_index, split_pos.clen, part.path_coords.back().clen);
1417 objects.push_back(obj2);
1418 }
1419 }
1420
1421 return objects;
1422 }
1423
changePathBounds(PathPartVector::size_type part_index,PathCoord::length_type start_len,PathCoord::length_type end_len)1424 void PathObject::changePathBounds(
1425 PathPartVector::size_type part_index,
1426 PathCoord::length_type start_len,
1427 PathCoord::length_type end_len)
1428 {
1429 update();
1430
1431 PathPart& part = path_parts[part_index];
1432 Q_ASSERT(end_len > start_len || part.isClosed());
1433
1434 auto part_size = part.size();
1435
1436 MapCoordVector out_coords;
1437 out_coords.reserve(part_size + 2);
1438
1439 auto part_begin = SplitPathCoord::begin(part.path_coords);
1440 auto part_end = SplitPathCoord::end(part.path_coords);
1441 if (start_len == part_end.clen && end_len == part_begin.clen)
1442 start_len = part_begin.clen;
1443 auto start = SplitPathCoord::at(start_len, part_begin);
1444
1445 if (end_len <= start_len)
1446 {
1447 if (start_len < part_end.clen)
1448 {
1449 // Make sure part_end has the right curve end points for start.
1450 part_end = SplitPathCoord::at(part_end.clen, start);
1451 part.copy(start, part_end, out_coords);
1452 out_coords.back().setHolePoint(false);
1453 out_coords.back().setClosePoint(false);
1454 }
1455
1456 if (end_len > part_begin.clen)
1457 {
1458 auto end = SplitPathCoord::at(end_len, part_begin);
1459 part.copy(part_begin, end, out_coords);
1460 }
1461
1462 Q_ASSERT(!out_coords.empty());
1463 }
1464 else
1465 {
1466 auto end = SplitPathCoord::at(end_len, start);
1467 part.copy(start, end, out_coords);
1468 }
1469
1470 out_coords.back().setHolePoint(true);
1471 out_coords.back().setClosePoint(false);
1472
1473 const auto copy_size = std::ptrdiff_t(std::min(out_coords.size(), std::size_t(part_size)));
1474 const auto part_start = begin(coords) + part.first_index;
1475 std::copy(begin(out_coords), begin(out_coords) + copy_size, part_start);
1476 if (copy_size < part_size)
1477 coords.erase(part_start + copy_size, part_start + part_size);
1478 else
1479 coords.insert(part_start + part_size, begin(out_coords) + copy_size, end(out_coords));
1480
1481 recalculateParts();
1482 setOutputDirty();
1483 }
1484
calcBezierPointDeletionRetainingShapeFactors(MapCoord p0,MapCoord p1,MapCoord p2,MapCoord q0,MapCoord q1,MapCoord q2,MapCoord q3,double & out_pfactor,double & out_qfactor)1485 void PathObject::calcBezierPointDeletionRetainingShapeFactors(MapCoord p0, MapCoord p1, MapCoord p2, MapCoord q0, MapCoord q1, MapCoord q2, MapCoord q3, double& out_pfactor, double& out_qfactor)
1486 {
1487 // Heuristic for the split parameter sp (zero to one)
1488 QBezier p_curve = QBezier::fromPoints(QPointF(p0), QPointF(p1), QPointF(p2), QPointF(q0));
1489 double p_length = p_curve.length(PathCoord::bezierError());
1490 QBezier q_curve = QBezier::fromPoints(QPointF(q0), QPointF(q1), QPointF(q2), QPointF(q3));
1491 double q_length = q_curve.length(PathCoord::bezierError());
1492 double sp = p_length / qMax(1e-08, p_length + q_length);
1493
1494 // Least squares curve fitting with the constraint that handles are on the same line as before.
1495 // To reproduce the formulas, run the following Matlab script:
1496 /*
1497
1498 clear all; close all; clc;
1499
1500 syms P0x P1x P2x P3x;
1501 syms Q1x Q2x Q3x;
1502 syms P0y P1y P2y P3y;
1503 syms Q1y Q2y Q3y;
1504
1505 syms u v w;
1506 syms pfactor qfactor;
1507
1508 Xx = (1-u)^3*P0x+3*(1-u)^2*u*P1x+3*(1-u)*u^2*P2x+u^3*P3x;
1509 Yx = (1-v)^3*P3x+3*(1-v)^2*v*Q1x+3*(1-v)*v^2*Q2x+v^3*Q3x;
1510 Zx = (1-w)^3*P0x+3*(1-w)^2*w*(P0x+pfactor*(P1x-P0x))+3*(1-w)*w^2*(Q3x+qfactor*(Q2x-Q3x))+w^3*Q3x;
1511
1512 Xy = (1-u)^3*P0y+3*(1-u)^2*u*P1y+3*(1-u)*u^2*P2y+u^3*P3y;
1513 Yy = (1-v)^3*P3y+3*(1-v)^2*v*Q1y+3*(1-v)*v^2*Q2y+v^3*Q3y;
1514 Zy = (1-w)^3*P0y+3*(1-w)^2*w*(P0y+pfactor*(P1y-P0y))+3*(1-w)*w^2*(Q3y+qfactor*(Q2y-Q3y))+w^3*Q3y;
1515
1516 syms sp; %split_param
1517 E = int((Zx - subs(Xx, u, w/sp))^2 + (Zy - subs(Xy, u, w/sp))^2, w, 0, sp) + int((Zx - subs(Yx, v, (w - sp) / (1 - sp)))^2 + (Zy - subs(Yy, v, (w - sp) / (1 - sp)))^2, w, sp, 1);
1518
1519 Epfactor = diff(E, pfactor);
1520 Eqfactor = diff(E, qfactor);
1521 S = solve(Epfactor, Eqfactor, pfactor, qfactor);
1522 S = [S.pfactor, S.qfactor]
1523
1524 */
1525
1526 double P0x = p0.x(), P1x = p1.x(), P2x = p2.x(), P3x = q0.x();
1527 double Q1x = q1.x(), Q2x = q2.x(), Q3x = q3.x();
1528 double P0y = p0.y(), P1y = p1.y(), P2y = p2.y(), P3y = q0.y();
1529 double Q1y = q1.y(), Q2y = q2.y(), Q3y = q3.y();
1530
1531 out_pfactor = -(126*P0x*pow(Q2x,3.0)*pow(sp,3.0) - 49*pow(P0x,2.0)*pow(Q3x,2.0) - 88*pow(P0x,2.0)*pow(Q2y,2.0) - 88*pow(P0y,2.0)*pow(Q2x,2.0) - 88*pow(P0x,2.0)*pow(Q3y,2.0) - 88*pow(P0y,2.0)*pow(Q3x,2.0) - 49*pow(P0y,2.0)*pow(Q2y,2.0) -
1532 49*pow(P0y,2.0)*pow(Q3y,2.0) - 49*pow(P0x,2.0)*pow(Q2x,2.0) + 21*P0x*pow(Q3x,3.0)*pow(sp,2.0) - 84*P0x*pow(Q2x,3.0)*pow(sp,4.0) + 14*P0x*pow(Q3x,3.0)*pow(sp,3.0) - 126*P1x*pow(Q2x,3.0)*pow(sp,3.0) -
1533 21*P1x*pow(Q3x,3.0)*pow(sp,2.0) - 21*P0x*pow(Q3x,3.0)*pow(sp,4.0) + 84*P1x*pow(Q2x,3.0)*pow(sp,4.0) - 14*P1x*pow(Q3x,3.0)*pow(sp,3.0) + 21*P1x*pow(Q3x,3.0)*pow(sp,4.0) + 126*P0y*pow(Q2y,3.0)*pow(sp,3.0) +
1534 21*P0y*pow(Q3y,3.0)*pow(sp,2.0) - 84*P0y*pow(Q2y,3.0)*pow(sp,4.0) + 14*P0y*pow(Q3y,3.0)*pow(sp,3.0) - 126*P1y*pow(Q2y,3.0)*pow(sp,3.0) - 21*P1y*pow(Q3y,3.0)*pow(sp,2.0) - 21*P0y*pow(Q3y,3.0)*pow(sp,4.0) +
1535 84*P1y*pow(Q2y,3.0)*pow(sp,4.0) - 14*P1y*pow(Q3y,3.0)*pow(sp,3.0) + 21*P1y*pow(Q3y,3.0)*pow(sp,4.0) + 84*pow(P0x,2.0)*pow(Q2x,2.0)*pow(sp,2.0) - 77*pow(P0x,2.0)*pow(Q2x,2.0)*pow(sp,3.0) +
1536 84*pow(P0x,2.0)*pow(Q3x,2.0)*pow(sp,2.0) - 168*pow(P1x,2.0)*pow(Q2x,2.0)*pow(sp,2.0) + 21*pow(P0x,2.0)*pow(Q2x,2.0)*pow(sp,4.0) - 77*pow(P0x,2.0)*pow(Q3x,2.0)*pow(sp,3.0) + 231*pow(P1x,2.0)*pow(Q2x,2.0)*pow(sp,3.0) -
1537 168*pow(P1x,2.0)*pow(Q3x,2.0)*pow(sp,2.0) + 21*pow(P0x,2.0)*pow(Q3x,2.0)*pow(sp,4.0) - 84*pow(P1x,2.0)*pow(Q2x,2.0)*pow(sp,4.0) + 231*pow(P1x,2.0)*pow(Q3x,2.0)*pow(sp,3.0) - 84*pow(P1x,2.0)*pow(Q3x,2.0)*pow(sp,4.0) +
1538 84*pow(P0x,2.0)*pow(Q2y,2.0)*pow(sp,2.0) + 84*pow(P0y,2.0)*pow(Q2x,2.0)*pow(sp,2.0) - 56*pow(P0x,2.0)*pow(Q2y,2.0)*pow(sp,3.0) + 84*pow(P0x,2.0)*pow(Q3y,2.0)*pow(sp,2.0) - 168*pow(P1x,2.0)*pow(Q2y,2.0)*pow(sp,2.0) -
1539 56*pow(P0y,2.0)*pow(Q2x,2.0)*pow(sp,3.0) + 84*pow(P0y,2.0)*pow(Q3x,2.0)*pow(sp,2.0) - 168*pow(P1y,2.0)*pow(Q2x,2.0)*pow(sp,2.0) + 12*pow(P0x,2.0)*pow(Q2y,2.0)*pow(sp,4.0) - 56*pow(P0x,2.0)*pow(Q3y,2.0)*pow(sp,3.0) +
1540 168*pow(P1x,2.0)*pow(Q2y,2.0)*pow(sp,3.0) - 168*pow(P1x,2.0)*pow(Q3y,2.0)*pow(sp,2.0) + 12*pow(P0y,2.0)*pow(Q2x,2.0)*pow(sp,4.0) - 56*pow(P0y,2.0)*pow(Q3x,2.0)*pow(sp,3.0) + 168*pow(P1y,2.0)*pow(Q2x,2.0)*pow(sp,3.0) -
1541 168*pow(P1y,2.0)*pow(Q3x,2.0)*pow(sp,2.0) + 12*pow(P0x,2.0)*pow(Q3y,2.0)*pow(sp,4.0) - 48*pow(P1x,2.0)*pow(Q2y,2.0)*pow(sp,4.0) + 168*pow(P1x,2.0)*pow(Q3y,2.0)*pow(sp,3.0) + 12*pow(P0y,2.0)*pow(Q3x,2.0)*pow(sp,4.0) -
1542 48*pow(P1y,2.0)*pow(Q2x,2.0)*pow(sp,4.0) + 168*pow(P1y,2.0)*pow(Q3x,2.0)*pow(sp,3.0) - 48*pow(P1x,2.0)*pow(Q3y,2.0)*pow(sp,4.0) - 48*pow(P1y,2.0)*pow(Q3x,2.0)*pow(sp,4.0) + 84*pow(P0y,2.0)*pow(Q2y,2.0)*pow(sp,2.0) -
1543 77*pow(P0y,2.0)*pow(Q2y,2.0)*pow(sp,3.0) + 84*pow(P0y,2.0)*pow(Q3y,2.0)*pow(sp,2.0) - 168*pow(P1y,2.0)*pow(Q2y,2.0)*pow(sp,2.0) + 21*pow(P0y,2.0)*pow(Q2y,2.0)*pow(sp,4.0) - 77*pow(P0y,2.0)*pow(Q3y,2.0)*pow(sp,3.0) +
1544 231*pow(P1y,2.0)*pow(Q2y,2.0)*pow(sp,3.0) - 168*pow(P1y,2.0)*pow(Q3y,2.0)*pow(sp,2.0) + 21*pow(P0y,2.0)*pow(Q3y,2.0)*pow(sp,4.0) - 84*pow(P1y,2.0)*pow(Q2y,2.0)*pow(sp,4.0) + 231*pow(P1y,2.0)*pow(Q3y,2.0)*pow(sp,3.0) -
1545 84*pow(P1y,2.0)*pow(Q3y,2.0)*pow(sp,4.0) + 49*P0x*P1x*pow(Q2x,2.0) + 49*P0x*P1x*pow(Q3x,2.0) + 28*P0x*P3x*pow(Q2x,2.0) + 28*P0x*P3x*pow(Q3x,2.0) - 28*P1x*P3x*pow(Q2x,2.0) - 28*P1x*P3x*pow(Q3x,2.0) + 88*P0x*P1x*pow(Q2y,2.0) +
1546 88*P0x*P1x*pow(Q3y,2.0) + 40*P0x*P3x*pow(Q2y,2.0) + 40*P0x*P3x*pow(Q3y,2.0) - 40*P1x*P3x*pow(Q2y,2.0) - 40*P1x*P3x*pow(Q3y,2.0) + 88*P0y*P1y*pow(Q2x,2.0) + 88*P0y*P1y*pow(Q3x,2.0) + 40*P0y*P3y*pow(Q2x,2.0) +
1547 40*P0y*P3y*pow(Q3x,2.0) - 40*P1y*P3y*pow(Q2x,2.0) - 40*P1y*P3y*pow(Q3x,2.0) + 49*P0y*P1y*pow(Q2y,2.0) + 49*P0y*P1y*pow(Q3y,2.0) + 28*P0y*P3y*pow(Q2y,2.0) + 28*P0y*P3y*pow(Q3y,2.0) - 28*P1y*P3y*pow(Q2y,2.0) -
1548 28*P1y*P3y*pow(Q3y,2.0) + 21*P0x*Q1x*pow(Q2x,2.0) + 21*P0x*Q1x*pow(Q3x,2.0) - 21*P1x*Q1x*pow(Q2x,2.0) - 21*P1x*Q1x*pow(Q3x,2.0) + 98*pow(P0x,2.0)*Q2x*Q3x + 48*P0x*Q1x*pow(Q2y,2.0) + 48*P0x*Q1x*pow(Q3y,2.0) -
1549 48*P1x*Q1x*pow(Q2y,2.0) - 48*P1x*Q1x*pow(Q3y,2.0) + 176*pow(P0y,2.0)*Q2x*Q3x + 48*P0y*pow(Q2x,2.0)*Q1y + 48*P0y*pow(Q3x,2.0)*Q1y - 48*P1y*pow(Q2x,2.0)*Q1y - 48*P1y*pow(Q3x,2.0)*Q1y + 176*pow(P0x,2.0)*Q2y*Q3y +
1550 21*P0y*Q1y*pow(Q2y,2.0) + 21*P0y*Q1y*pow(Q3y,2.0) - 21*P1y*Q1y*pow(Q2y,2.0) - 21*P1y*Q1y*pow(Q3y,2.0) + 98*pow(P0y,2.0)*Q2y*Q3y - 42*P0x*pow(Q2x,3.0)*sp + 42*P1x*pow(Q2x,3.0)*sp - 42*P0y*pow(Q2y,3.0)*sp +
1551 42*P1y*pow(Q2y,3.0)*sp + 84*P0x*P3x*pow(Q2x,2.0)*sp + 84*P0x*P3x*pow(Q3x,2.0)*sp - 84*P1x*P3x*pow(Q2x,2.0)*sp - 84*P1x*P3x*pow(Q3x,2.0)*sp + 120*P0x*P3x*pow(Q2y,2.0)*sp + 120*P0x*P3x*pow(Q3y,2.0)*sp -
1552 120*P1x*P3x*pow(Q2y,2.0)*sp - 120*P1x*P3x*pow(Q3y,2.0)*sp + 120*P0y*P3y*pow(Q2x,2.0)*sp + 120*P0y*P3y*pow(Q3x,2.0)*sp - 120*P1y*P3y*pow(Q2x,2.0)*sp - 120*P1y*P3y*pow(Q3x,2.0)*sp + 84*P0y*P3y*pow(Q2y,2.0)*sp +
1553 84*P0y*P3y*pow(Q3y,2.0)*sp - 84*P1y*P3y*pow(Q2y,2.0)*sp - 84*P1y*P3y*pow(Q3y,2.0)*sp - 42*P0x*Q1x*pow(Q2x,2.0)*sp - 42*P0x*Q1x*pow(Q3x,2.0)*sp + 42*P1x*Q1x*pow(Q2x,2.0)*sp - 42*P0x*Q2x*pow(Q3x,2.0)*sp +
1554 84*P0x*pow(Q2x,2.0)*Q3x*sp + 42*P1x*Q1x*pow(Q3x,2.0)*sp + 42*P1x*Q2x*pow(Q3x,2.0)*sp - 84*P1x*pow(Q2x,2.0)*Q3x*sp - 24*P0x*Q1x*pow(Q2y,2.0)*sp - 24*P0x*Q1x*pow(Q3y,2.0)*sp - 42*P0x*Q2x*pow(Q2y,2.0)*sp +
1555 24*P1x*Q1x*pow(Q2y,2.0)*sp - 96*P0x*Q2x*pow(Q3y,2.0)*sp - 54*P0x*Q3x*pow(Q2y,2.0)*sp + 24*P1x*Q1x*pow(Q3y,2.0)*sp + 42*P1x*Q2x*pow(Q2y,2.0)*sp + 96*P1x*Q2x*pow(Q3y,2.0)*sp + 54*P1x*Q3x*pow(Q2y,2.0)*sp -
1556 24*P0y*pow(Q2x,2.0)*Q1y*sp - 42*P0y*pow(Q2x,2.0)*Q2y*sp - 24*P0y*pow(Q3x,2.0)*Q1y*sp + 24*P1y*pow(Q2x,2.0)*Q1y*sp - 54*P0y*pow(Q2x,2.0)*Q3y*sp - 96*P0y*pow(Q3x,2.0)*Q2y*sp + 42*P1y*pow(Q2x,2.0)*Q2y*sp +
1557 24*P1y*pow(Q3x,2.0)*Q1y*sp + 54*P1y*pow(Q2x,2.0)*Q3y*sp + 96*P1y*pow(Q3x,2.0)*Q2y*sp - 42*P0y*Q1y*pow(Q2y,2.0)*sp - 42*P0y*Q1y*pow(Q3y,2.0)*sp + 42*P1y*Q1y*pow(Q2y,2.0)*sp - 42*P0y*Q2y*pow(Q3y,2.0)*sp +
1558 84*P0y*pow(Q2y,2.0)*Q3y*sp + 42*P1y*Q1y*pow(Q3y,2.0)*sp + 42*P1y*Q2y*pow(Q3y,2.0)*sp - 84*P1y*pow(Q2y,2.0)*Q3y*sp + 84*P0x*P1x*pow(Q2x,2.0)*pow(sp,2.0) - 154*P0x*P1x*pow(Q2x,2.0)*pow(sp,3.0) +
1559 84*P0x*P1x*pow(Q3x,2.0)*pow(sp,2.0) + 252*P0x*P2x*pow(Q2x,2.0)*pow(sp,2.0) + 63*P0x*P1x*pow(Q2x,2.0)*pow(sp,4.0) - 154*P0x*P1x*pow(Q3x,2.0)*pow(sp,3.0) - 462*P0x*P2x*pow(Q2x,2.0)*pow(sp,3.0) +
1560 252*P0x*P2x*pow(Q3x,2.0)*pow(sp,2.0) - 336*P0x*P3x*pow(Q2x,2.0)*pow(sp,2.0) - 252*P1x*P2x*pow(Q2x,2.0)*pow(sp,2.0) + 63*P0x*P1x*pow(Q3x,2.0)*pow(sp,4.0) + 210*P0x*P2x*pow(Q2x,2.0)*pow(sp,4.0) -
1561 462*P0x*P2x*pow(Q3x,2.0)*pow(sp,3.0) + 210*P0x*P3x*pow(Q2x,2.0)*pow(sp,3.0) - 336*P0x*P3x*pow(Q3x,2.0)*pow(sp,2.0) + 462*P1x*P2x*pow(Q2x,2.0)*pow(sp,3.0) - 252*P1x*P2x*pow(Q3x,2.0)*pow(sp,2.0) +
1562 336*P1x*P3x*pow(Q2x,2.0)*pow(sp,2.0) + 210*P0x*P2x*pow(Q3x,2.0)*pow(sp,4.0) + 210*P0x*P3x*pow(Q3x,2.0)*pow(sp,3.0) - 210*P1x*P2x*pow(Q2x,2.0)*pow(sp,4.0) + 462*P1x*P2x*pow(Q3x,2.0)*pow(sp,3.0) -
1563 210*P1x*P3x*pow(Q2x,2.0)*pow(sp,3.0) + 336*P1x*P3x*pow(Q3x,2.0)*pow(sp,2.0) - 210*P1x*P2x*pow(Q3x,2.0)*pow(sp,4.0) - 210*P1x*P3x*pow(Q3x,2.0)*pow(sp,3.0) + 84*P0x*P1x*pow(Q2y,2.0)*pow(sp,2.0) -
1564 112*P0x*P1x*pow(Q2y,2.0)*pow(sp,3.0) + 84*P0x*P1x*pow(Q3y,2.0)*pow(sp,2.0) + 252*P0x*P2x*pow(Q2y,2.0)*pow(sp,2.0) + 36*P0x*P1x*pow(Q2y,2.0)*pow(sp,4.0) - 112*P0x*P1x*pow(Q3y,2.0)*pow(sp,3.0) -
1565 336*P0x*P2x*pow(Q2y,2.0)*pow(sp,3.0) + 252*P0x*P2x*pow(Q3y,2.0)*pow(sp,2.0) - 264*P0x*P3x*pow(Q2y,2.0)*pow(sp,2.0) - 252*P1x*P2x*pow(Q2y,2.0)*pow(sp,2.0) + 36*P0x*P1x*pow(Q3y,2.0)*pow(sp,4.0) +
1566 120*P0x*P2x*pow(Q2y,2.0)*pow(sp,4.0) - 336*P0x*P2x*pow(Q3y,2.0)*pow(sp,3.0) + 120*P0x*P3x*pow(Q2y,2.0)*pow(sp,3.0) - 264*P0x*P3x*pow(Q3y,2.0)*pow(sp,2.0) + 336*P1x*P2x*pow(Q2y,2.0)*pow(sp,3.0) -
1567 252*P1x*P2x*pow(Q3y,2.0)*pow(sp,2.0) + 264*P1x*P3x*pow(Q2y,2.0)*pow(sp,2.0) + 120*P0x*P2x*pow(Q3y,2.0)*pow(sp,4.0) + 120*P0x*P3x*pow(Q3y,2.0)*pow(sp,3.0) - 120*P1x*P2x*pow(Q2y,2.0)*pow(sp,4.0) +
1568 336*P1x*P2x*pow(Q3y,2.0)*pow(sp,3.0) - 120*P1x*P3x*pow(Q2y,2.0)*pow(sp,3.0) + 264*P1x*P3x*pow(Q3y,2.0)*pow(sp,2.0) - 120*P1x*P2x*pow(Q3y,2.0)*pow(sp,4.0) - 120*P1x*P3x*pow(Q3y,2.0)*pow(sp,3.0) +
1569 84*P0y*P1y*pow(Q2x,2.0)*pow(sp,2.0) - 112*P0y*P1y*pow(Q2x,2.0)*pow(sp,3.0) + 84*P0y*P1y*pow(Q3x,2.0)*pow(sp,2.0) + 252*P0y*P2y*pow(Q2x,2.0)*pow(sp,2.0) + 36*P0y*P1y*pow(Q2x,2.0)*pow(sp,4.0) -
1570 112*P0y*P1y*pow(Q3x,2.0)*pow(sp,3.0) - 336*P0y*P2y*pow(Q2x,2.0)*pow(sp,3.0) + 252*P0y*P2y*pow(Q3x,2.0)*pow(sp,2.0) - 264*P0y*P3y*pow(Q2x,2.0)*pow(sp,2.0) - 252*P1y*P2y*pow(Q2x,2.0)*pow(sp,2.0) +
1571 36*P0y*P1y*pow(Q3x,2.0)*pow(sp,4.0) + 120*P0y*P2y*pow(Q2x,2.0)*pow(sp,4.0) - 336*P0y*P2y*pow(Q3x,2.0)*pow(sp,3.0) + 120*P0y*P3y*pow(Q2x,2.0)*pow(sp,3.0) - 264*P0y*P3y*pow(Q3x,2.0)*pow(sp,2.0) +
1572 336*P1y*P2y*pow(Q2x,2.0)*pow(sp,3.0) - 252*P1y*P2y*pow(Q3x,2.0)*pow(sp,2.0) + 264*P1y*P3y*pow(Q2x,2.0)*pow(sp,2.0) + 120*P0y*P2y*pow(Q3x,2.0)*pow(sp,4.0) + 120*P0y*P3y*pow(Q3x,2.0)*pow(sp,3.0) -
1573 120*P1y*P2y*pow(Q2x,2.0)*pow(sp,4.0) + 336*P1y*P2y*pow(Q3x,2.0)*pow(sp,3.0) - 120*P1y*P3y*pow(Q2x,2.0)*pow(sp,3.0) + 264*P1y*P3y*pow(Q3x,2.0)*pow(sp,2.0) - 120*P1y*P2y*pow(Q3x,2.0)*pow(sp,4.0) -
1574 120*P1y*P3y*pow(Q3x,2.0)*pow(sp,3.0) + 84*P0y*P1y*pow(Q2y,2.0)*pow(sp,2.0) - 154*P0y*P1y*pow(Q2y,2.0)*pow(sp,3.0) + 84*P0y*P1y*pow(Q3y,2.0)*pow(sp,2.0) + 252*P0y*P2y*pow(Q2y,2.0)*pow(sp,2.0) +
1575 63*P0y*P1y*pow(Q2y,2.0)*pow(sp,4.0) - 154*P0y*P1y*pow(Q3y,2.0)*pow(sp,3.0) - 462*P0y*P2y*pow(Q2y,2.0)*pow(sp,3.0) + 252*P0y*P2y*pow(Q3y,2.0)*pow(sp,2.0) - 336*P0y*P3y*pow(Q2y,2.0)*pow(sp,2.0) -
1576 252*P1y*P2y*pow(Q2y,2.0)*pow(sp,2.0) + 63*P0y*P1y*pow(Q3y,2.0)*pow(sp,4.0) + 210*P0y*P2y*pow(Q2y,2.0)*pow(sp,4.0) - 462*P0y*P2y*pow(Q3y,2.0)*pow(sp,3.0) + 210*P0y*P3y*pow(Q2y,2.0)*pow(sp,3.0) -
1577 336*P0y*P3y*pow(Q3y,2.0)*pow(sp,2.0) + 462*P1y*P2y*pow(Q2y,2.0)*pow(sp,3.0) - 252*P1y*P2y*pow(Q3y,2.0)*pow(sp,2.0) + 336*P1y*P3y*pow(Q2y,2.0)*pow(sp,2.0) + 210*P0y*P2y*pow(Q3y,2.0)*pow(sp,4.0) +
1578 210*P0y*P3y*pow(Q3y,2.0)*pow(sp,3.0) - 210*P1y*P2y*pow(Q2y,2.0)*pow(sp,4.0) + 462*P1y*P2y*pow(Q3y,2.0)*pow(sp,3.0) - 210*P1y*P3y*pow(Q2y,2.0)*pow(sp,3.0) + 336*P1y*P3y*pow(Q3y,2.0)*pow(sp,2.0) -
1579 210*P1y*P2y*pow(Q3y,2.0)*pow(sp,4.0) - 210*P1y*P3y*pow(Q3y,2.0)*pow(sp,3.0) - 189*P0x*Q1x*pow(Q2x,2.0)*pow(sp,2.0) + 420*P0x*Q1x*pow(Q2x,2.0)*pow(sp,3.0) - 189*P0x*Q1x*pow(Q3x,2.0)*pow(sp,2.0) +
1580 189*P1x*Q1x*pow(Q2x,2.0)*pow(sp,2.0) - 210*P0x*Q1x*pow(Q2x,2.0)*pow(sp,4.0) + 420*P0x*Q1x*pow(Q3x,2.0)*pow(sp,3.0) - 42*P0x*Q2x*pow(Q3x,2.0)*pow(sp,2.0) + 21*P0x*pow(Q2x,2.0)*Q3x*pow(sp,2.0) -
1581 420*P1x*Q1x*pow(Q2x,2.0)*pow(sp,3.0) + 189*P1x*Q1x*pow(Q3x,2.0)*pow(sp,2.0) - 168*pow(P0x,2.0)*Q2x*Q3x*pow(sp,2.0) - 210*P0x*Q1x*pow(Q3x,2.0)*pow(sp,4.0) + 98*P0x*Q2x*pow(Q3x,2.0)*pow(sp,3.0) -
1582 238*P0x*pow(Q2x,2.0)*Q3x*pow(sp,3.0) + 210*P1x*Q1x*pow(Q2x,2.0)*pow(sp,4.0) - 420*P1x*Q1x*pow(Q3x,2.0)*pow(sp,3.0) + 42*P1x*Q2x*pow(Q3x,2.0)*pow(sp,2.0) - 21*P1x*pow(Q2x,2.0)*Q3x*pow(sp,2.0) +
1583 154*pow(P0x,2.0)*Q2x*Q3x*pow(sp,3.0) + 336*pow(P1x,2.0)*Q2x*Q3x*pow(sp,2.0) - 42*P0x*Q2x*pow(Q3x,2.0)*pow(sp,4.0) + 147*P0x*pow(Q2x,2.0)*Q3x*pow(sp,4.0) + 210*P1x*Q1x*pow(Q3x,2.0)*pow(sp,4.0) -
1584 98*P1x*Q2x*pow(Q3x,2.0)*pow(sp,3.0) + 238*P1x*pow(Q2x,2.0)*Q3x*pow(sp,3.0) - 42*pow(P0x,2.0)*Q2x*Q3x*pow(sp,4.0) - 462*pow(P1x,2.0)*Q2x*Q3x*pow(sp,3.0) + 42*P1x*Q2x*pow(Q3x,2.0)*pow(sp,4.0) -
1585 147*P1x*pow(Q2x,2.0)*Q3x*pow(sp,4.0) + 168*pow(P1x,2.0)*Q2x*Q3x*pow(sp,4.0) - 216*P0x*Q1x*pow(Q2y,2.0)*pow(sp,2.0) + 312*P0x*Q1x*pow(Q2y,2.0)*pow(sp,3.0) - 216*P0x*Q1x*pow(Q3y,2.0)*pow(sp,2.0) +
1586 216*P1x*Q1x*pow(Q2y,2.0)*pow(sp,2.0) - 120*P0x*Q1x*pow(Q2y,2.0)*pow(sp,4.0) + 312*P0x*Q1x*pow(Q3y,2.0)*pow(sp,3.0) + 126*P0x*Q2x*pow(Q2y,2.0)*pow(sp,3.0) - 45*P0x*Q2x*pow(Q3y,2.0)*pow(sp,2.0) -
1587 24*P0x*Q3x*pow(Q2y,2.0)*pow(sp,2.0) - 312*P1x*Q1x*pow(Q2y,2.0)*pow(sp,3.0) + 216*P1x*Q1x*pow(Q3y,2.0)*pow(sp,2.0) - 168*pow(P0y,2.0)*Q2x*Q3x*pow(sp,2.0) - 120*P0x*Q1x*pow(Q3y,2.0)*pow(sp,4.0) -
1588 84*P0x*Q2x*pow(Q2y,2.0)*pow(sp,4.0) + 114*P0x*Q2x*pow(Q3y,2.0)*pow(sp,3.0) + 2*P0x*Q3x*pow(Q2y,2.0)*pow(sp,3.0) + 21*P0x*Q3x*pow(Q3y,2.0)*pow(sp,2.0) + 120*P1x*Q1x*pow(Q2y,2.0)*pow(sp,4.0) -
1589 312*P1x*Q1x*pow(Q3y,2.0)*pow(sp,3.0) - 126*P1x*Q2x*pow(Q2y,2.0)*pow(sp,3.0) + 45*P1x*Q2x*pow(Q3y,2.0)*pow(sp,2.0) + 24*P1x*Q3x*pow(Q2y,2.0)*pow(sp,2.0) + 112*pow(P0y,2.0)*Q2x*Q3x*pow(sp,3.0) +
1590 336*pow(P1y,2.0)*Q2x*Q3x*pow(sp,2.0) - 39*P0x*Q2x*pow(Q3y,2.0)*pow(sp,4.0) + 24*P0x*Q3x*pow(Q2y,2.0)*pow(sp,4.0) + 14*P0x*Q3x*pow(Q3y,2.0)*pow(sp,3.0) + 120*P1x*Q1x*pow(Q3y,2.0)*pow(sp,4.0) +
1591 84*P1x*Q2x*pow(Q2y,2.0)*pow(sp,4.0) - 114*P1x*Q2x*pow(Q3y,2.0)*pow(sp,3.0) - 2*P1x*Q3x*pow(Q2y,2.0)*pow(sp,3.0) - 21*P1x*Q3x*pow(Q3y,2.0)*pow(sp,2.0) - 24*pow(P0y,2.0)*Q2x*Q3x*pow(sp,4.0) -
1592 336*pow(P1y,2.0)*Q2x*Q3x*pow(sp,3.0) - 21*P0x*Q3x*pow(Q3y,2.0)*pow(sp,4.0) + 39*P1x*Q2x*pow(Q3y,2.0)*pow(sp,4.0) - 24*P1x*Q3x*pow(Q2y,2.0)*pow(sp,4.0) - 14*P1x*Q3x*pow(Q3y,2.0)*pow(sp,3.0) +
1593 96*pow(P1y,2.0)*Q2x*Q3x*pow(sp,4.0) + 21*P1x*Q3x*pow(Q3y,2.0)*pow(sp,4.0) - 216*P0y*pow(Q2x,2.0)*Q1y*pow(sp,2.0) + 312*P0y*pow(Q2x,2.0)*Q1y*pow(sp,3.0) - 216*P0y*pow(Q3x,2.0)*Q1y*pow(sp,2.0) +
1594 216*P1y*pow(Q2x,2.0)*Q1y*pow(sp,2.0) - 120*P0y*pow(Q2x,2.0)*Q1y*pow(sp,4.0) + 126*P0y*pow(Q2x,2.0)*Q2y*pow(sp,3.0) - 24*P0y*pow(Q2x,2.0)*Q3y*pow(sp,2.0) + 312*P0y*pow(Q3x,2.0)*Q1y*pow(sp,3.0) -
1595 45*P0y*pow(Q3x,2.0)*Q2y*pow(sp,2.0) - 312*P1y*pow(Q2x,2.0)*Q1y*pow(sp,3.0) + 216*P1y*pow(Q3x,2.0)*Q1y*pow(sp,2.0) - 168*pow(P0x,2.0)*Q2y*Q3y*pow(sp,2.0) - 84*P0y*pow(Q2x,2.0)*Q2y*pow(sp,4.0) +
1596 2*P0y*pow(Q2x,2.0)*Q3y*pow(sp,3.0) - 120*P0y*pow(Q3x,2.0)*Q1y*pow(sp,4.0) + 114*P0y*pow(Q3x,2.0)*Q2y*pow(sp,3.0) + 21*P0y*pow(Q3x,2.0)*Q3y*pow(sp,2.0) + 120*P1y*pow(Q2x,2.0)*Q1y*pow(sp,4.0) -
1597 126*P1y*pow(Q2x,2.0)*Q2y*pow(sp,3.0) + 24*P1y*pow(Q2x,2.0)*Q3y*pow(sp,2.0) - 312*P1y*pow(Q3x,2.0)*Q1y*pow(sp,3.0) + 45*P1y*pow(Q3x,2.0)*Q2y*pow(sp,2.0) + 112*pow(P0x,2.0)*Q2y*Q3y*pow(sp,3.0) +
1598 336*pow(P1x,2.0)*Q2y*Q3y*pow(sp,2.0) + 24*P0y*pow(Q2x,2.0)*Q3y*pow(sp,4.0) - 39*P0y*pow(Q3x,2.0)*Q2y*pow(sp,4.0) + 14*P0y*pow(Q3x,2.0)*Q3y*pow(sp,3.0) + 84*P1y*pow(Q2x,2.0)*Q2y*pow(sp,4.0) -
1599 2*P1y*pow(Q2x,2.0)*Q3y*pow(sp,3.0) + 120*P1y*pow(Q3x,2.0)*Q1y*pow(sp,4.0) - 114*P1y*pow(Q3x,2.0)*Q2y*pow(sp,3.0) - 21*P1y*pow(Q3x,2.0)*Q3y*pow(sp,2.0) - 24*pow(P0x,2.0)*Q2y*Q3y*pow(sp,4.0) -
1600 336*pow(P1x,2.0)*Q2y*Q3y*pow(sp,3.0) - 21*P0y*pow(Q3x,2.0)*Q3y*pow(sp,4.0) - 24*P1y*pow(Q2x,2.0)*Q3y*pow(sp,4.0) + 39*P1y*pow(Q3x,2.0)*Q2y*pow(sp,4.0) - 14*P1y*pow(Q3x,2.0)*Q3y*pow(sp,3.0) +
1601 96*pow(P1x,2.0)*Q2y*Q3y*pow(sp,4.0) + 21*P1y*pow(Q3x,2.0)*Q3y*pow(sp,4.0) - 189*P0y*Q1y*pow(Q2y,2.0)*pow(sp,2.0) + 420*P0y*Q1y*pow(Q2y,2.0)*pow(sp,3.0) - 189*P0y*Q1y*pow(Q3y,2.0)*pow(sp,2.0) +
1602 189*P1y*Q1y*pow(Q2y,2.0)*pow(sp,2.0) - 210*P0y*Q1y*pow(Q2y,2.0)*pow(sp,4.0) + 420*P0y*Q1y*pow(Q3y,2.0)*pow(sp,3.0) - 42*P0y*Q2y*pow(Q3y,2.0)*pow(sp,2.0) + 21*P0y*pow(Q2y,2.0)*Q3y*pow(sp,2.0) -
1603 420*P1y*Q1y*pow(Q2y,2.0)*pow(sp,3.0) + 189*P1y*Q1y*pow(Q3y,2.0)*pow(sp,2.0) - 168*pow(P0y,2.0)*Q2y*Q3y*pow(sp,2.0) - 210*P0y*Q1y*pow(Q3y,2.0)*pow(sp,4.0) + 98*P0y*Q2y*pow(Q3y,2.0)*pow(sp,3.0) -
1604 238*P0y*pow(Q2y,2.0)*Q3y*pow(sp,3.0) + 210*P1y*Q1y*pow(Q2y,2.0)*pow(sp,4.0) - 420*P1y*Q1y*pow(Q3y,2.0)*pow(sp,3.0) + 42*P1y*Q2y*pow(Q3y,2.0)*pow(sp,2.0) - 21*P1y*pow(Q2y,2.0)*Q3y*pow(sp,2.0) +
1605 154*pow(P0y,2.0)*Q2y*Q3y*pow(sp,3.0) + 336*pow(P1y,2.0)*Q2y*Q3y*pow(sp,2.0) - 42*P0y*Q2y*pow(Q3y,2.0)*pow(sp,4.0) + 147*P0y*pow(Q2y,2.0)*Q3y*pow(sp,4.0) + 210*P1y*Q1y*pow(Q3y,2.0)*pow(sp,4.0) -
1606 98*P1y*Q2y*pow(Q3y,2.0)*pow(sp,3.0) + 238*P1y*pow(Q2y,2.0)*Q3y*pow(sp,3.0) - 42*pow(P0y,2.0)*Q2y*Q3y*pow(sp,4.0) - 462*pow(P1y,2.0)*Q2y*Q3y*pow(sp,3.0) + 42*P1y*Q2y*pow(Q3y,2.0)*pow(sp,4.0) -
1607 147*P1y*pow(Q2y,2.0)*Q3y*pow(sp,4.0) + 168*pow(P1y,2.0)*Q2y*Q3y*pow(sp,4.0) - 98*P0x*P1x*Q2x*Q3x - 56*P0x*P3x*Q2x*Q3x + 56*P1x*P3x*Q2x*Q3x + 78*P0x*P0y*Q2x*Q2y - 78*P0x*P0y*Q2x*Q3y - 78*P0x*P0y*Q3x*Q2y -
1608 39*P0x*P1y*Q2x*Q2y - 39*P1x*P0y*Q2x*Q2y - 176*P0x*P1x*Q2y*Q3y + 78*P0x*P0y*Q3x*Q3y + 39*P0x*P1y*Q2x*Q3y + 39*P0x*P1y*Q3x*Q2y + 39*P1x*P0y*Q2x*Q3y + 39*P1x*P0y*Q3x*Q2y - 176*P0y*P1y*Q2x*Q3x -
1609 39*P0x*P1y*Q3x*Q3y - 12*P0x*P3y*Q2x*Q2y - 39*P1x*P0y*Q3x*Q3y - 12*P3x*P0y*Q2x*Q2y - 80*P0x*P3x*Q2y*Q3y + 12*P0x*P3y*Q2x*Q3y + 12*P0x*P3y*Q3x*Q2y + 12*P1x*P3y*Q2x*Q2y + 12*P3x*P0y*Q2x*Q3y +
1610 12*P3x*P0y*Q3x*Q2y + 12*P3x*P1y*Q2x*Q2y - 80*P0y*P3y*Q2x*Q3x - 12*P0x*P3y*Q3x*Q3y + 80*P1x*P3x*Q2y*Q3y - 12*P1x*P3y*Q2x*Q3y - 12*P1x*P3y*Q3x*Q2y - 12*P3x*P0y*Q3x*Q3y - 12*P3x*P1y*Q2x*Q3y -
1611 12*P3x*P1y*Q3x*Q2y + 80*P1y*P3y*Q2x*Q3x + 12*P1x*P3y*Q3x*Q3y + 12*P3x*P1y*Q3x*Q3y - 98*P0y*P1y*Q2y*Q3y - 56*P0y*P3y*Q2y*Q3y + 56*P1y*P3y*Q2y*Q3y - 42*P0x*Q1x*Q2x*Q3x + 42*P1x*Q1x*Q2x*Q3x - 27*P0x*Q2x*Q1y*Q2y -
1612 27*P0y*Q1x*Q2x*Q2y - 96*P0x*Q1x*Q2y*Q3y + 27*P0x*Q2x*Q1y*Q3y + 27*P0x*Q3x*Q1y*Q2y + 27*P1x*Q2x*Q1y*Q2y + 27*P0y*Q1x*Q2x*Q3y + 27*P0y*Q1x*Q3x*Q2y - 96*P0y*Q2x*Q3x*Q1y + 27*P1y*Q1x*Q2x*Q2y - 27*P0x*Q3x*Q1y*Q3y +
1613 96*P1x*Q1x*Q2y*Q3y - 27*P1x*Q2x*Q1y*Q3y - 27*P1x*Q3x*Q1y*Q2y - 27*P0y*Q1x*Q3x*Q3y - 27*P1y*Q1x*Q2x*Q3y - 27*P1y*Q1x*Q3x*Q2y + 96*P1y*Q2x*Q3x*Q1y + 27*P1x*Q3x*Q1y*Q3y + 27*P1y*Q1x*Q3x*Q3y - 42*P0y*Q1y*Q2y*Q3y +
1614 42*P1y*Q1y*Q2y*Q3y - 168*P0x*P3x*Q2x*Q3x*sp + 168*P1x*P3x*Q2x*Q3x*sp - 36*P0x*P3y*Q2x*Q2y*sp - 36*P3x*P0y*Q2x*Q2y*sp - 240*P0x*P3x*Q2y*Q3y*sp + 36*P0x*P3y*Q2x*Q3y*sp + 36*P0x*P3y*Q3x*Q2y*sp + 36*P1x*P3y*Q2x*Q2y*sp +
1615 36*P3x*P0y*Q2x*Q3y*sp + 36*P3x*P0y*Q3x*Q2y*sp + 36*P3x*P1y*Q2x*Q2y*sp - 240*P0y*P3y*Q2x*Q3x*sp - 36*P0x*P3y*Q3x*Q3y*sp + 240*P1x*P3x*Q2y*Q3y*sp - 36*P1x*P3y*Q2x*Q3y*sp - 36*P1x*P3y*Q3x*Q2y*sp - 36*P3x*P0y*Q3x*Q3y*sp -
1616 36*P3x*P1y*Q2x*Q3y*sp - 36*P3x*P1y*Q3x*Q2y*sp + 240*P1y*P3y*Q2x*Q3x*sp + 36*P1x*P3y*Q3x*Q3y*sp + 36*P3x*P1y*Q3x*Q3y*sp - 168*P0y*P3y*Q2y*Q3y*sp + 168*P1y*P3y*Q2y*Q3y*sp + 84*P0x*Q1x*Q2x*Q3x*sp - 84*P1x*Q1x*Q2x*Q3x*sp -
1617 18*P0x*Q2x*Q1y*Q2y*sp - 18*P0y*Q1x*Q2x*Q2y*sp + 48*P0x*Q1x*Q2y*Q3y*sp + 18*P0x*Q2x*Q1y*Q3y*sp + 18*P0x*Q3x*Q1y*Q2y*sp + 18*P1x*Q2x*Q1y*Q2y*sp + 18*P0y*Q1x*Q2x*Q3y*sp + 18*P0y*Q1x*Q3x*Q2y*sp + 48*P0y*Q2x*Q3x*Q1y*sp +
1618 18*P1y*Q1x*Q2x*Q2y*sp + 138*P0x*Q2x*Q2y*Q3y*sp - 18*P0x*Q3x*Q1y*Q3y*sp - 48*P1x*Q1x*Q2y*Q3y*sp - 18*P1x*Q2x*Q1y*Q3y*sp - 18*P1x*Q3x*Q1y*Q2y*sp - 18*P0y*Q1x*Q3x*Q3y*sp + 138*P0y*Q2x*Q3x*Q2y*sp - 18*P1y*Q1x*Q2x*Q3y*sp -
1619 18*P1y*Q1x*Q3x*Q2y*sp - 48*P1y*Q2x*Q3x*Q1y*sp + 54*P0x*Q3x*Q2y*Q3y*sp - 138*P1x*Q2x*Q2y*Q3y*sp + 18*P1x*Q3x*Q1y*Q3y*sp + 54*P0y*Q2x*Q3x*Q3y*sp + 18*P1y*Q1x*Q3x*Q3y*sp - 138*P1y*Q2x*Q3x*Q2y*sp - 54*P1x*Q3x*Q2y*Q3y*sp -
1620 54*P1y*Q2x*Q3x*Q3y*sp + 84*P0y*Q1y*Q2y*Q3y*sp - 84*P1y*Q1y*Q2y*Q3y*sp - 168*P0x*P1x*Q2x*Q3x*pow(sp,2.0) + 308*P0x*P1x*Q2x*Q3x*pow(sp,3.0) - 504*P0x*P2x*Q2x*Q3x*pow(sp,2.0) - 126*P0x*P1x*Q2x*Q3x*pow(sp,4.0) +
1621 924*P0x*P2x*Q2x*Q3x*pow(sp,3.0) + 672*P0x*P3x*Q2x*Q3x*pow(sp,2.0) + 504*P1x*P2x*Q2x*Q3x*pow(sp,2.0) - 420*P0x*P2x*Q2x*Q3x*pow(sp,4.0) - 420*P0x*P3x*Q2x*Q3x*pow(sp,3.0) - 924*P1x*P2x*Q2x*Q3x*pow(sp,3.0) -
1622 672*P1x*P3x*Q2x*Q3x*pow(sp,2.0) + 420*P1x*P2x*Q2x*Q3x*pow(sp,4.0) + 420*P1x*P3x*Q2x*Q3x*pow(sp,3.0) - 42*P0x*P0y*Q2x*Q2y*pow(sp,3.0) - 168*P0x*P1x*Q2y*Q3y*pow(sp,2.0) + 18*P0x*P0y*Q2x*Q2y*pow(sp,4.0) +
1623 42*P0x*P0y*Q2x*Q3y*pow(sp,3.0) + 42*P0x*P0y*Q3x*Q2y*pow(sp,3.0) - 42*P0x*P1y*Q2x*Q2y*pow(sp,3.0) - 42*P1x*P0y*Q2x*Q2y*pow(sp,3.0) - 168*P0y*P1y*Q2x*Q3x*pow(sp,2.0) + 224*P0x*P1x*Q2y*Q3y*pow(sp,3.0) -
1624 504*P0x*P2x*Q2y*Q3y*pow(sp,2.0) - 18*P0x*P0y*Q2x*Q3y*pow(sp,4.0) - 18*P0x*P0y*Q3x*Q2y*pow(sp,4.0) - 42*P0x*P0y*Q3x*Q3y*pow(sp,3.0) + 27*P0x*P1y*Q2x*Q2y*pow(sp,4.0) + 42*P0x*P1y*Q2x*Q3y*pow(sp,3.0) +
1625 42*P0x*P1y*Q3x*Q2y*pow(sp,3.0) - 126*P0x*P2y*Q2x*Q2y*pow(sp,3.0) - 72*P0x*P3y*Q2x*Q2y*pow(sp,2.0) + 27*P1x*P0y*Q2x*Q2y*pow(sp,4.0) + 42*P1x*P0y*Q2x*Q3y*pow(sp,3.0) + 42*P1x*P0y*Q3x*Q2y*pow(sp,3.0) +
1626 126*P1x*P1y*Q2x*Q2y*pow(sp,3.0) - 126*P2x*P0y*Q2x*Q2y*pow(sp,3.0) - 72*P3x*P0y*Q2x*Q2y*pow(sp,2.0) + 224*P0y*P1y*Q2x*Q3x*pow(sp,3.0) - 504*P0y*P2y*Q2x*Q3x*pow(sp,2.0) - 72*P0x*P1x*Q2y*Q3y*pow(sp,4.0) +
1627 672*P0x*P2x*Q2y*Q3y*pow(sp,3.0) + 528*P0x*P3x*Q2y*Q3y*pow(sp,2.0) + 18*P0x*P0y*Q3x*Q3y*pow(sp,4.0) - 27*P0x*P1y*Q2x*Q3y*pow(sp,4.0) - 27*P0x*P1y*Q3x*Q2y*pow(sp,4.0) - 42*P0x*P1y*Q3x*Q3y*pow(sp,3.0) +
1628 90*P0x*P2y*Q2x*Q2y*pow(sp,4.0) + 126*P0x*P2y*Q2x*Q3y*pow(sp,3.0) + 126*P0x*P2y*Q3x*Q2y*pow(sp,3.0) + 90*P0x*P3y*Q2x*Q2y*pow(sp,3.0) + 72*P0x*P3y*Q2x*Q3y*pow(sp,2.0) + 72*P0x*P3y*Q3x*Q2y*pow(sp,2.0) +
1629 504*P1x*P2x*Q2y*Q3y*pow(sp,2.0) - 27*P1x*P0y*Q2x*Q3y*pow(sp,4.0) - 27*P1x*P0y*Q3x*Q2y*pow(sp,4.0) - 42*P1x*P0y*Q3x*Q3y*pow(sp,3.0) - 72*P1x*P1y*Q2x*Q2y*pow(sp,4.0) - 126*P1x*P1y*Q2x*Q3y*pow(sp,3.0) -
1630 126*P1x*P1y*Q3x*Q2y*pow(sp,3.0) + 126*P1x*P2y*Q2x*Q2y*pow(sp,3.0) + 72*P1x*P3y*Q2x*Q2y*pow(sp,2.0) + 90*P2x*P0y*Q2x*Q2y*pow(sp,4.0) + 126*P2x*P0y*Q2x*Q3y*pow(sp,3.0) + 126*P2x*P0y*Q3x*Q2y*pow(sp,3.0) +
1631 126*P2x*P1y*Q2x*Q2y*pow(sp,3.0) + 90*P3x*P0y*Q2x*Q2y*pow(sp,3.0) + 72*P3x*P0y*Q2x*Q3y*pow(sp,2.0) + 72*P3x*P0y*Q3x*Q2y*pow(sp,2.0) + 72*P3x*P1y*Q2x*Q2y*pow(sp,2.0) - 72*P0y*P1y*Q2x*Q3x*pow(sp,4.0) +
1632 672*P0y*P2y*Q2x*Q3x*pow(sp,3.0) + 528*P0y*P3y*Q2x*Q3x*pow(sp,2.0) + 504*P1y*P2y*Q2x*Q3x*pow(sp,2.0) - 240*P0x*P2x*Q2y*Q3y*pow(sp,4.0) - 240*P0x*P3x*Q2y*Q3y*pow(sp,3.0) + 27*P0x*P1y*Q3x*Q3y*pow(sp,4.0) -
1633 90*P0x*P2y*Q2x*Q3y*pow(sp,4.0) - 90*P0x*P2y*Q3x*Q2y*pow(sp,4.0) - 126*P0x*P2y*Q3x*Q3y*pow(sp,3.0) - 90*P0x*P3y*Q2x*Q3y*pow(sp,3.0) - 90*P0x*P3y*Q3x*Q2y*pow(sp,3.0) - 72*P0x*P3y*Q3x*Q3y*pow(sp,2.0) -
1634 672*P1x*P2x*Q2y*Q3y*pow(sp,3.0) - 528*P1x*P3x*Q2y*Q3y*pow(sp,2.0) + 27*P1x*P0y*Q3x*Q3y*pow(sp,4.0) + 72*P1x*P1y*Q2x*Q3y*pow(sp,4.0) + 72*P1x*P1y*Q3x*Q2y*pow(sp,4.0) + 126*P1x*P1y*Q3x*Q3y*pow(sp,3.0) -
1635 90*P1x*P2y*Q2x*Q2y*pow(sp,4.0) - 126*P1x*P2y*Q2x*Q3y*pow(sp,3.0) - 126*P1x*P2y*Q3x*Q2y*pow(sp,3.0) - 90*P1x*P3y*Q2x*Q2y*pow(sp,3.0) - 72*P1x*P3y*Q2x*Q3y*pow(sp,2.0) - 72*P1x*P3y*Q3x*Q2y*pow(sp,2.0) -
1636 90*P2x*P0y*Q2x*Q3y*pow(sp,4.0) - 90*P2x*P0y*Q3x*Q2y*pow(sp,4.0) - 126*P2x*P0y*Q3x*Q3y*pow(sp,3.0) - 90*P2x*P1y*Q2x*Q2y*pow(sp,4.0) - 126*P2x*P1y*Q2x*Q3y*pow(sp,3.0) - 126*P2x*P1y*Q3x*Q2y*pow(sp,3.0) -
1637 90*P3x*P0y*Q2x*Q3y*pow(sp,3.0) - 90*P3x*P0y*Q3x*Q2y*pow(sp,3.0) - 72*P3x*P0y*Q3x*Q3y*pow(sp,2.0) - 90*P3x*P1y*Q2x*Q2y*pow(sp,3.0) - 72*P3x*P1y*Q2x*Q3y*pow(sp,2.0) - 72*P3x*P1y*Q3x*Q2y*pow(sp,2.0) -
1638 240*P0y*P2y*Q2x*Q3x*pow(sp,4.0) - 240*P0y*P3y*Q2x*Q3x*pow(sp,3.0) - 672*P1y*P2y*Q2x*Q3x*pow(sp,3.0) - 528*P1y*P3y*Q2x*Q3x*pow(sp,2.0) + 90*P0x*P2y*Q3x*Q3y*pow(sp,4.0) + 90*P0x*P3y*Q3x*Q3y*pow(sp,3.0) +
1639 240*P1x*P2x*Q2y*Q3y*pow(sp,4.0) + 240*P1x*P3x*Q2y*Q3y*pow(sp,3.0) - 72*P1x*P1y*Q3x*Q3y*pow(sp,4.0) + 90*P1x*P2y*Q2x*Q3y*pow(sp,4.0) + 90*P1x*P2y*Q3x*Q2y*pow(sp,4.0) + 126*P1x*P2y*Q3x*Q3y*pow(sp,3.0) +
1640 90*P1x*P3y*Q2x*Q3y*pow(sp,3.0) + 90*P1x*P3y*Q3x*Q2y*pow(sp,3.0) + 72*P1x*P3y*Q3x*Q3y*pow(sp,2.0) + 90*P2x*P0y*Q3x*Q3y*pow(sp,4.0) + 90*P2x*P1y*Q2x*Q3y*pow(sp,4.0) + 90*P2x*P1y*Q3x*Q2y*pow(sp,4.0) +
1641 126*P2x*P1y*Q3x*Q3y*pow(sp,3.0) + 90*P3x*P0y*Q3x*Q3y*pow(sp,3.0) + 90*P3x*P1y*Q2x*Q3y*pow(sp,3.0) + 90*P3x*P1y*Q3x*Q2y*pow(sp,3.0) + 72*P3x*P1y*Q3x*Q3y*pow(sp,2.0) + 240*P1y*P2y*Q2x*Q3x*pow(sp,4.0) +
1642 240*P1y*P3y*Q2x*Q3x*pow(sp,3.0) - 90*P1x*P2y*Q3x*Q3y*pow(sp,4.0) - 90*P1x*P3y*Q3x*Q3y*pow(sp,3.0) - 90*P2x*P1y*Q3x*Q3y*pow(sp,4.0) - 90*P3x*P1y*Q3x*Q3y*pow(sp,3.0) - 168*P0y*P1y*Q2y*Q3y*pow(sp,2.0) +
1643 308*P0y*P1y*Q2y*Q3y*pow(sp,3.0) - 504*P0y*P2y*Q2y*Q3y*pow(sp,2.0) - 126*P0y*P1y*Q2y*Q3y*pow(sp,4.0) + 924*P0y*P2y*Q2y*Q3y*pow(sp,3.0) + 672*P0y*P3y*Q2y*Q3y*pow(sp,2.0) + 504*P1y*P2y*Q2y*Q3y*pow(sp,2.0) -
1644 420*P0y*P2y*Q2y*Q3y*pow(sp,4.0) - 420*P0y*P3y*Q2y*Q3y*pow(sp,3.0) - 924*P1y*P2y*Q2y*Q3y*pow(sp,3.0) - 672*P1y*P3y*Q2y*Q3y*pow(sp,2.0) + 420*P1y*P2y*Q2y*Q3y*pow(sp,4.0) + 420*P1y*P3y*Q2y*Q3y*pow(sp,3.0) +
1645 378*P0x*Q1x*Q2x*Q3x*pow(sp,2.0) - 840*P0x*Q1x*Q2x*Q3x*pow(sp,3.0) - 378*P1x*Q1x*Q2x*Q3x*pow(sp,2.0) + 420*P0x*Q1x*Q2x*Q3x*pow(sp,4.0) + 840*P1x*Q1x*Q2x*Q3x*pow(sp,3.0) - 420*P1x*Q1x*Q2x*Q3x*pow(sp,4.0) +
1646 27*P0x*Q2x*Q1y*Q2y*pow(sp,2.0) + 27*P0y*Q1x*Q2x*Q2y*pow(sp,2.0) + 432*P0x*Q1x*Q2y*Q3y*pow(sp,2.0) + 108*P0x*Q2x*Q1y*Q2y*pow(sp,3.0) - 27*P0x*Q2x*Q1y*Q3y*pow(sp,2.0) - 27*P0x*Q3x*Q1y*Q2y*pow(sp,2.0) -
1647 27*P1x*Q2x*Q1y*Q2y*pow(sp,2.0) + 108*P0y*Q1x*Q2x*Q2y*pow(sp,3.0) - 27*P0y*Q1x*Q2x*Q3y*pow(sp,2.0) - 27*P0y*Q1x*Q3x*Q2y*pow(sp,2.0) + 432*P0y*Q2x*Q3x*Q1y*pow(sp,2.0) - 27*P1y*Q1x*Q2x*Q2y*pow(sp,2.0) -
1648 624*P0x*Q1x*Q2y*Q3y*pow(sp,3.0) - 90*P0x*Q2x*Q1y*Q2y*pow(sp,4.0) - 108*P0x*Q2x*Q1y*Q3y*pow(sp,3.0) + 45*P0x*Q2x*Q2y*Q3y*pow(sp,2.0) - 108*P0x*Q3x*Q1y*Q2y*pow(sp,3.0) + 27*P0x*Q3x*Q1y*Q3y*pow(sp,2.0) -
1649 432*P1x*Q1x*Q2y*Q3y*pow(sp,2.0) - 108*P1x*Q2x*Q1y*Q2y*pow(sp,3.0) + 27*P1x*Q2x*Q1y*Q3y*pow(sp,2.0) + 27*P1x*Q3x*Q1y*Q2y*pow(sp,2.0) - 90*P0y*Q1x*Q2x*Q2y*pow(sp,4.0) - 108*P0y*Q1x*Q2x*Q3y*pow(sp,3.0) -
1650 108*P0y*Q1x*Q3x*Q2y*pow(sp,3.0) + 27*P0y*Q1x*Q3x*Q3y*pow(sp,2.0) - 624*P0y*Q2x*Q3x*Q1y*pow(sp,3.0) + 45*P0y*Q2x*Q3x*Q2y*pow(sp,2.0) - 108*P1y*Q1x*Q2x*Q2y*pow(sp,3.0) + 27*P1y*Q1x*Q2x*Q3y*pow(sp,2.0) +
1651 27*P1y*Q1x*Q3x*Q2y*pow(sp,2.0) - 432*P1y*Q2x*Q3x*Q1y*pow(sp,2.0) + 240*P0x*Q1x*Q2y*Q3y*pow(sp,4.0) + 90*P0x*Q2x*Q1y*Q3y*pow(sp,4.0) - 240*P0x*Q2x*Q2y*Q3y*pow(sp,3.0) + 90*P0x*Q3x*Q1y*Q2y*pow(sp,4.0) +
1652 108*P0x*Q3x*Q1y*Q3y*pow(sp,3.0) + 3*P0x*Q3x*Q2y*Q3y*pow(sp,2.0) + 624*P1x*Q1x*Q2y*Q3y*pow(sp,3.0) + 90*P1x*Q2x*Q1y*Q2y*pow(sp,4.0) + 108*P1x*Q2x*Q1y*Q3y*pow(sp,3.0) - 45*P1x*Q2x*Q2y*Q3y*pow(sp,2.0) +
1653 108*P1x*Q3x*Q1y*Q2y*pow(sp,3.0) - 27*P1x*Q3x*Q1y*Q3y*pow(sp,2.0) + 90*P0y*Q1x*Q2x*Q3y*pow(sp,4.0) + 90*P0y*Q1x*Q3x*Q2y*pow(sp,4.0) + 108*P0y*Q1x*Q3x*Q3y*pow(sp,3.0) + 240*P0y*Q2x*Q3x*Q1y*pow(sp,4.0) -
1654 240*P0y*Q2x*Q3x*Q2y*pow(sp,3.0) + 3*P0y*Q2x*Q3x*Q3y*pow(sp,2.0) + 90*P1y*Q1x*Q2x*Q2y*pow(sp,4.0) + 108*P1y*Q1x*Q2x*Q3y*pow(sp,3.0) + 108*P1y*Q1x*Q3x*Q2y*pow(sp,3.0) - 27*P1y*Q1x*Q3x*Q3y*pow(sp,2.0) +
1655 624*P1y*Q2x*Q3x*Q1y*pow(sp,3.0) - 45*P1y*Q2x*Q3x*Q2y*pow(sp,2.0) + 123*P0x*Q2x*Q2y*Q3y*pow(sp,4.0) - 90*P0x*Q3x*Q1y*Q3y*pow(sp,4.0) - 16*P0x*Q3x*Q2y*Q3y*pow(sp,3.0) - 240*P1x*Q1x*Q2y*Q3y*pow(sp,4.0) -
1656 90*P1x*Q2x*Q1y*Q3y*pow(sp,4.0) + 240*P1x*Q2x*Q2y*Q3y*pow(sp,3.0) - 90*P1x*Q3x*Q1y*Q2y*pow(sp,4.0) - 108*P1x*Q3x*Q1y*Q3y*pow(sp,3.0) - 3*P1x*Q3x*Q2y*Q3y*pow(sp,2.0) - 90*P0y*Q1x*Q3x*Q3y*pow(sp,4.0) +
1657 123*P0y*Q2x*Q3x*Q2y*pow(sp,4.0) - 16*P0y*Q2x*Q3x*Q3y*pow(sp,3.0) - 90*P1y*Q1x*Q2x*Q3y*pow(sp,4.0) - 90*P1y*Q1x*Q3x*Q2y*pow(sp,4.0) - 108*P1y*Q1x*Q3x*Q3y*pow(sp,3.0) - 240*P1y*Q2x*Q3x*Q1y*pow(sp,4.0) +
1658 240*P1y*Q2x*Q3x*Q2y*pow(sp,3.0) - 3*P1y*Q2x*Q3x*Q3y*pow(sp,2.0) - 3*P0x*Q3x*Q2y*Q3y*pow(sp,4.0) - 123*P1x*Q2x*Q2y*Q3y*pow(sp,4.0) + 90*P1x*Q3x*Q1y*Q3y*pow(sp,4.0) + 16*P1x*Q3x*Q2y*Q3y*pow(sp,3.0) -
1659 3*P0y*Q2x*Q3x*Q3y*pow(sp,4.0) + 90*P1y*Q1x*Q3x*Q3y*pow(sp,4.0) - 123*P1y*Q2x*Q3x*Q2y*pow(sp,4.0) + 16*P1y*Q2x*Q3x*Q3y*pow(sp,3.0) + 3*P1x*Q3x*Q2y*Q3y*pow(sp,4.0) + 3*P1y*Q2x*Q3x*Q3y*pow(sp,4.0) +
1660 378*P0y*Q1y*Q2y*Q3y*pow(sp,2.0) - 840*P0y*Q1y*Q2y*Q3y*pow(sp,3.0) - 378*P1y*Q1y*Q2y*Q3y*pow(sp,2.0) + 420*P0y*Q1y*Q2y*Q3y*pow(sp,4.0) + 840*P1y*Q1y*Q2y*Q3y*pow(sp,3.0) -
1661 420*P1y*Q1y*Q2y*Q3y*pow(sp,4.0))/(3*(7*pow(P0x,2.0)*pow(Q2x,2.0) + 7*pow(P0x,2.0)*pow(Q3x,2.0) + 7*pow(P1x,2.0)*pow(Q2x,2.0) + 7*pow(P1x,2.0)*pow(Q3x,2.0) + 16*pow(P0x,2.0)*pow(Q2y,2.0) +
1662 16*pow(P0y,2.0)*pow(Q2x,2.0) + 16*pow(P0x,2.0)*pow(Q3y,2.0) + 16*pow(P1x,2.0)*pow(Q2y,2.0) + 16*pow(P0y,2.0)*pow(Q3x,2.0) + 16*pow(P1y,2.0)*pow(Q2x,2.0) + 16*pow(P1x,2.0)*pow(Q3y,2.0) +
1663 16*pow(P1y,2.0)*pow(Q3x,2.0) + 7*pow(P0y,2.0)*pow(Q2y,2.0) + 7*pow(P0y,2.0)*pow(Q3y,2.0) + 7*pow(P1y,2.0)*pow(Q2y,2.0) + 7*pow(P1y,2.0)*pow(Q3y,2.0) - 14*P0x*P1x*pow(Q2x,2.0) - 14*P0x*P1x*pow(Q3x,2.0) -
1664 32*P0x*P1x*pow(Q2y,2.0) - 32*P0x*P1x*pow(Q3y,2.0) - 32*P0y*P1y*pow(Q2x,2.0) - 32*P0y*P1y*pow(Q3x,2.0) - 14*P0y*P1y*pow(Q2y,2.0) - 14*P0y*P1y*pow(Q3y,2.0) - 14*pow(P0x,2.0)*Q2x*Q3x - 14*pow(P1x,2.0)*Q2x*Q3x -
1665 32*pow(P0y,2.0)*Q2x*Q3x - 32*pow(P1y,2.0)*Q2x*Q3x - 32*pow(P0x,2.0)*Q2y*Q3y - 32*pow(P1x,2.0)*Q2y*Q3y - 14*pow(P0y,2.0)*Q2y*Q3y - 14*pow(P1y,2.0)*Q2y*Q3y + 28*P0x*P1x*Q2x*Q3x - 18*P0x*P0y*Q2x*Q2y + 18*P0x*P0y*Q2x*Q3y +
1666 18*P0x*P0y*Q3x*Q2y + 18*P0x*P1y*Q2x*Q2y + 18*P1x*P0y*Q2x*Q2y + 64*P0x*P1x*Q2y*Q3y - 18*P0x*P0y*Q3x*Q3y - 18*P0x*P1y*Q2x*Q3y - 18*P0x*P1y*Q3x*Q2y - 18*P1x*P0y*Q2x*Q3y - 18*P1x*P0y*Q3x*Q2y - 18*P1x*P1y*Q2x*Q2y +
1667 64*P0y*P1y*Q2x*Q3x + 18*P0x*P1y*Q3x*Q3y + 18*P1x*P0y*Q3x*Q3y + 18*P1x*P1y*Q2x*Q3y + 18*P1x*P1y*Q3x*Q2y - 18*P1x*P1y*Q3x*Q3y + 28*P0y*P1y*Q2y*Q3y));
1668
1669 out_qfactor = (14*pow(P0x,3.0)*Q2x - 14*pow(P0x,3.0)*Q3x + 14*pow(P0y,3.0)*Q2y - 14*pow(P0y,3.0)*Q3y + 21*pow(P0x,2.0)*pow(Q2x,2.0) + 21*pow(P0x,2.0)*pow(Q3x,2.0) + 21*pow(P1x,2.0)*pow(Q2x,2.0) + 21*pow(P1x,2.0)*pow(Q3x,2.0) +
1670 48*pow(P0x,2.0)*pow(Q2y,2.0) + 48*pow(P0y,2.0)*pow(Q2x,2.0) + 48*pow(P0x,2.0)*pow(Q3y,2.0) + 48*pow(P1x,2.0)*pow(Q2y,2.0) + 48*pow(P0y,2.0)*pow(Q3x,2.0) + 48*pow(P1y,2.0)*pow(Q2x,2.0) + 48*pow(P1x,2.0)*pow(Q3y,2.0) +
1671 48*pow(P1y,2.0)*pow(Q3x,2.0) + 21*pow(P0y,2.0)*pow(Q2y,2.0) + 21*pow(P0y,2.0)*pow(Q3y,2.0) + 21*pow(P1y,2.0)*pow(Q2y,2.0) + 21*pow(P1y,2.0)*pow(Q3y,2.0) + 21*pow(P0x,2.0)*pow(Q2x,2.0)*sp + 21*pow(P0x,2.0)*pow(Q3x,2.0)*sp -
1672 63*pow(P0x,3.0)*Q2x*pow(sp,2.0) + 21*pow(P1x,2.0)*pow(Q2x,2.0)*sp + 70*pow(P0x,3.0)*Q2x*pow(sp,3.0) + 63*pow(P0x,3.0)*Q3x*pow(sp,2.0) + 21*pow(P1x,2.0)*pow(Q3x,2.0)*sp - 126*pow(P1x,3.0)*Q2x*pow(sp,2.0) -
1673 21*pow(P0x,3.0)*Q2x*pow(sp,4.0) - 70*pow(P0x,3.0)*Q3x*pow(sp,3.0) + 210*pow(P1x,3.0)*Q2x*pow(sp,3.0) + 126*pow(P1x,3.0)*Q3x*pow(sp,2.0) + 21*pow(P0x,3.0)*Q3x*pow(sp,4.0) - 84*pow(P1x,3.0)*Q2x*pow(sp,4.0) -
1674 210*pow(P1x,3.0)*Q3x*pow(sp,3.0) + 84*pow(P1x,3.0)*Q3x*pow(sp,4.0) - 24*pow(P0x,2.0)*pow(Q2y,2.0)*sp - 24*pow(P0y,2.0)*pow(Q2x,2.0)*sp + 48*pow(P0x,2.0)*pow(Q3y,2.0)*sp - 24*pow(P1x,2.0)*pow(Q2y,2.0)*sp +
1675 48*pow(P0y,2.0)*pow(Q3x,2.0)*sp - 24*pow(P1y,2.0)*pow(Q2x,2.0)*sp + 48*pow(P1x,2.0)*pow(Q3y,2.0)*sp + 48*pow(P1y,2.0)*pow(Q3x,2.0)*sp + 21*pow(P0y,2.0)*pow(Q2y,2.0)*sp + 21*pow(P0y,2.0)*pow(Q3y,2.0)*sp -
1676 63*pow(P0y,3.0)*Q2y*pow(sp,2.0) + 21*pow(P1y,2.0)*pow(Q2y,2.0)*sp + 70*pow(P0y,3.0)*Q2y*pow(sp,3.0) + 63*pow(P0y,3.0)*Q3y*pow(sp,2.0) + 21*pow(P1y,2.0)*pow(Q3y,2.0)*sp - 126*pow(P1y,3.0)*Q2y*pow(sp,2.0) -
1677 21*pow(P0y,3.0)*Q2y*pow(sp,4.0) - 70*pow(P0y,3.0)*Q3y*pow(sp,3.0) + 210*pow(P1y,3.0)*Q2y*pow(sp,3.0) + 126*pow(P1y,3.0)*Q3y*pow(sp,2.0) + 21*pow(P0y,3.0)*Q3y*pow(sp,4.0) - 84*pow(P1y,3.0)*Q2y*pow(sp,4.0) -
1678 210*pow(P1y,3.0)*Q3y*pow(sp,3.0) + 84*pow(P1y,3.0)*Q3y*pow(sp,4.0) - 21*pow(P0x,2.0)*pow(Q2x,2.0)*pow(sp,2.0) - 105*pow(P0x,2.0)*pow(Q2x,2.0)*pow(sp,3.0) + 21*pow(P0x,2.0)*pow(Q3x,2.0)*pow(sp,2.0) -
1679 21*pow(P1x,2.0)*pow(Q2x,2.0)*pow(sp,2.0) + 84*pow(P0x,2.0)*pow(Q2x,2.0)*pow(sp,4.0) + 7*pow(P0x,2.0)*pow(Q3x,2.0)*pow(sp,3.0) - 105*pow(P1x,2.0)*pow(Q2x,2.0)*pow(sp,3.0) + 21*pow(P1x,2.0)*pow(Q3x,2.0)*pow(sp,2.0) -
1680 21*pow(P0x,2.0)*pow(Q3x,2.0)*pow(sp,4.0) + 84*pow(P1x,2.0)*pow(Q2x,2.0)*pow(sp,4.0) + 7*pow(P1x,2.0)*pow(Q3x,2.0)*pow(sp,3.0) - 21*pow(P1x,2.0)*pow(Q3x,2.0)*pow(sp,4.0) - 48*pow(P0x,2.0)*pow(Q2y,2.0)*pow(sp,2.0) -
1681 48*pow(P0y,2.0)*pow(Q2x,2.0)*pow(sp,2.0) - 24*pow(P0x,2.0)*pow(Q2y,2.0)*pow(sp,3.0) + 12*pow(P0x,2.0)*pow(Q3y,2.0)*pow(sp,2.0) - 48*pow(P1x,2.0)*pow(Q2y,2.0)*pow(sp,2.0) - 24*pow(P0y,2.0)*pow(Q2x,2.0)*pow(sp,3.0) +
1682 12*pow(P0y,2.0)*pow(Q3x,2.0)*pow(sp,2.0) - 48*pow(P1y,2.0)*pow(Q2x,2.0)*pow(sp,2.0) + 48*pow(P0x,2.0)*pow(Q2y,2.0)*pow(sp,4.0) - 8*pow(P0x,2.0)*pow(Q3y,2.0)*pow(sp,3.0) - 24*pow(P1x,2.0)*pow(Q2y,2.0)*pow(sp,3.0) +
1683 12*pow(P1x,2.0)*pow(Q3y,2.0)*pow(sp,2.0) + 48*pow(P0y,2.0)*pow(Q2x,2.0)*pow(sp,4.0) - 8*pow(P0y,2.0)*pow(Q3x,2.0)*pow(sp,3.0) - 24*pow(P1y,2.0)*pow(Q2x,2.0)*pow(sp,3.0) + 12*pow(P1y,2.0)*pow(Q3x,2.0)*pow(sp,2.0) -
1684 12*pow(P0x,2.0)*pow(Q3y,2.0)*pow(sp,4.0) + 48*pow(P1x,2.0)*pow(Q2y,2.0)*pow(sp,4.0) - 8*pow(P1x,2.0)*pow(Q3y,2.0)*pow(sp,3.0) - 12*pow(P0y,2.0)*pow(Q3x,2.0)*pow(sp,4.0) + 48*pow(P1y,2.0)*pow(Q2x,2.0)*pow(sp,4.0) -
1685 8*pow(P1y,2.0)*pow(Q3x,2.0)*pow(sp,3.0) - 12*pow(P1x,2.0)*pow(Q3y,2.0)*pow(sp,4.0) - 12*pow(P1y,2.0)*pow(Q3x,2.0)*pow(sp,4.0) - 21*pow(P0y,2.0)*pow(Q2y,2.0)*pow(sp,2.0) - 105*pow(P0y,2.0)*pow(Q2y,2.0)*pow(sp,3.0) +
1686 21*pow(P0y,2.0)*pow(Q3y,2.0)*pow(sp,2.0) - 21*pow(P1y,2.0)*pow(Q2y,2.0)*pow(sp,2.0) + 84*pow(P0y,2.0)*pow(Q2y,2.0)*pow(sp,4.0) + 7*pow(P0y,2.0)*pow(Q3y,2.0)*pow(sp,3.0) - 105*pow(P1y,2.0)*pow(Q2y,2.0)*pow(sp,3.0) +
1687 21*pow(P1y,2.0)*pow(Q3y,2.0)*pow(sp,2.0) - 21*pow(P0y,2.0)*pow(Q3y,2.0)*pow(sp,4.0) + 84*pow(P1y,2.0)*pow(Q2y,2.0)*pow(sp,4.0) + 7*pow(P1y,2.0)*pow(Q3y,2.0)*pow(sp,3.0) - 21*pow(P1y,2.0)*pow(Q3y,2.0)*pow(sp,4.0) -
1688 42*P0x*P1x*pow(Q2x,2.0) + 14*P0x*pow(P1x,2.0)*Q2x - 28*pow(P0x,2.0)*P1x*Q2x - 42*P0x*P1x*pow(Q3x,2.0) - 14*P0x*pow(P1x,2.0)*Q3x + 28*pow(P0x,2.0)*P1x*Q3x - 14*pow(P0x,2.0)*P3x*Q2x + 14*pow(P0x,2.0)*P3x*Q3x -
1689 14*pow(P1x,2.0)*P3x*Q2x + 14*pow(P1x,2.0)*P3x*Q3x + 14*P0x*pow(P0y,2.0)*Q2x - 96*P0x*P1x*pow(Q2y,2.0) - 14*P0x*pow(P0y,2.0)*Q3x - 52*P0x*pow(P1y,2.0)*Q2x - 66*P1x*pow(P0y,2.0)*Q2x - 96*P0x*P1x*pow(Q3y,2.0) +
1690 52*P0x*pow(P1y,2.0)*Q3x + 66*P1x*pow(P0y,2.0)*Q3x + 16*P3x*pow(P0y,2.0)*Q2x - 16*P3x*pow(P0y,2.0)*Q3x + 16*P3x*pow(P1y,2.0)*Q2x - 16*P3x*pow(P1y,2.0)*Q3x + 14*pow(P0x,2.0)*P0y*Q2y - 96*P0y*P1y*pow(Q2x,2.0) -
1691 14*pow(P0x,2.0)*P0y*Q3y - 66*pow(P0x,2.0)*P1y*Q2y - 52*pow(P1x,2.0)*P0y*Q2y - 96*P0y*P1y*pow(Q3x,2.0) + 66*pow(P0x,2.0)*P1y*Q3y + 52*pow(P1x,2.0)*P0y*Q3y + 16*pow(P0x,2.0)*P3y*Q2y - 16*pow(P0x,2.0)*P3y*Q3y +
1692 16*pow(P1x,2.0)*P3y*Q2y - 16*pow(P1x,2.0)*P3y*Q3y - 42*P0y*P1y*pow(Q2y,2.0) + 14*P0y*pow(P1y,2.0)*Q2y - 28*pow(P0y,2.0)*P1y*Q2y - 42*P0y*P1y*pow(Q3y,2.0) - 14*P0y*pow(P1y,2.0)*Q3y + 28*pow(P0y,2.0)*P1y*Q3y -
1693 14*pow(P0y,2.0)*P3y*Q2y + 14*pow(P0y,2.0)*P3y*Q3y - 14*pow(P1y,2.0)*P3y*Q2y + 14*pow(P1y,2.0)*P3y*Q3y - 42*pow(P0x,2.0)*Q2x*Q3x - 42*pow(P1x,2.0)*Q2x*Q3x + 36*pow(P0y,2.0)*Q1x*Q2x - 36*pow(P0y,2.0)*Q1x*Q3x +
1694 36*pow(P1y,2.0)*Q1x*Q2x - 96*pow(P0y,2.0)*Q2x*Q3x - 36*pow(P1y,2.0)*Q1x*Q3x - 96*pow(P1y,2.0)*Q2x*Q3x + 36*pow(P0x,2.0)*Q1y*Q2y - 36*pow(P0x,2.0)*Q1y*Q3y + 36*pow(P1x,2.0)*Q1y*Q2y - 96*pow(P0x,2.0)*Q2y*Q3y -
1695 36*pow(P1x,2.0)*Q1y*Q3y - 96*pow(P1x,2.0)*Q2y*Q3y - 42*pow(P0y,2.0)*Q2y*Q3y - 42*pow(P1y,2.0)*Q2y*Q3y - 42*P0x*P1x*pow(Q2x,2.0)*sp - 42*P0x*P1x*pow(Q3x,2.0)*sp - 42*pow(P0x,2.0)*P3x*Q2x*sp + 42*pow(P0x,2.0)*P3x*Q3x*sp -
1696 42*pow(P1x,2.0)*P3x*Q2x*sp + 42*pow(P1x,2.0)*P3x*Q3x*sp + 48*P0x*P1x*pow(Q2y,2.0)*sp - 96*P0x*P1x*pow(Q3y,2.0)*sp + 48*P3x*pow(P0y,2.0)*Q2x*sp - 48*P3x*pow(P0y,2.0)*Q3x*sp + 48*P3x*pow(P1y,2.0)*Q2x*sp -
1697 48*P3x*pow(P1y,2.0)*Q3x*sp + 48*P0y*P1y*pow(Q2x,2.0)*sp - 96*P0y*P1y*pow(Q3x,2.0)*sp + 48*pow(P0x,2.0)*P3y*Q2y*sp - 48*pow(P0x,2.0)*P3y*Q3y*sp + 48*pow(P1x,2.0)*P3y*Q2y*sp - 48*pow(P1x,2.0)*P3y*Q3y*sp -
1698 42*P0y*P1y*pow(Q2y,2.0)*sp - 42*P0y*P1y*pow(Q3y,2.0)*sp - 42*pow(P0y,2.0)*P3y*Q2y*sp + 42*pow(P0y,2.0)*P3y*Q3y*sp - 42*pow(P1y,2.0)*P3y*Q2y*sp + 42*pow(P1y,2.0)*P3y*Q3y*sp + 42*pow(P0x,2.0)*Q1x*Q2x*sp -
1699 42*pow(P0x,2.0)*Q1x*Q3x*sp + 42*pow(P1x,2.0)*Q1x*Q2x*sp - 42*pow(P0x,2.0)*Q2x*Q3x*sp - 42*pow(P1x,2.0)*Q1x*Q3x*sp - 42*pow(P1x,2.0)*Q2x*Q3x*sp + 24*pow(P0y,2.0)*Q1x*Q2x*sp - 24*pow(P0y,2.0)*Q1x*Q3x*sp +
1700 24*pow(P1y,2.0)*Q1x*Q2x*sp - 24*pow(P0y,2.0)*Q2x*Q3x*sp - 24*pow(P1y,2.0)*Q1x*Q3x*sp - 24*pow(P1y,2.0)*Q2x*Q3x*sp + 24*pow(P0x,2.0)*Q1y*Q2y*sp - 24*pow(P0x,2.0)*Q1y*Q3y*sp + 24*pow(P1x,2.0)*Q1y*Q2y*sp -
1701 24*pow(P0x,2.0)*Q2y*Q3y*sp - 24*pow(P1x,2.0)*Q1y*Q3y*sp - 24*pow(P1x,2.0)*Q2y*Q3y*sp + 42*pow(P0y,2.0)*Q1y*Q2y*sp - 42*pow(P0y,2.0)*Q1y*Q3y*sp + 42*pow(P1y,2.0)*Q1y*Q2y*sp - 42*pow(P0y,2.0)*Q2y*Q3y*sp -
1702 42*pow(P1y,2.0)*Q1y*Q3y*sp - 42*pow(P1y,2.0)*Q2y*Q3y*sp + 42*P0x*P1x*pow(Q2x,2.0)*pow(sp,2.0) + 189*P0x*pow(P1x,2.0)*Q2x*pow(sp,2.0) + 210*P0x*P1x*pow(Q2x,2.0)*pow(sp,3.0) - 42*P0x*P1x*pow(Q3x,2.0)*pow(sp,2.0) -
1703 350*P0x*pow(P1x,2.0)*Q2x*pow(sp,3.0) - 189*P0x*pow(P1x,2.0)*Q3x*pow(sp,2.0) + 70*pow(P0x,2.0)*P1x*Q2x*pow(sp,3.0) - 189*pow(P0x,2.0)*P2x*Q2x*pow(sp,2.0) - 168*P0x*P1x*pow(Q2x,2.0)*pow(sp,4.0) -
1704 14*P0x*P1x*pow(Q3x,2.0)*pow(sp,3.0) + 147*P0x*pow(P1x,2.0)*Q2x*pow(sp,4.0) + 350*P0x*pow(P1x,2.0)*Q3x*pow(sp,3.0) - 42*pow(P0x,2.0)*P1x*Q2x*pow(sp,4.0) - 70*pow(P0x,2.0)*P1x*Q3x*pow(sp,3.0) +
1705 420*pow(P0x,2.0)*P2x*Q2x*pow(sp,3.0) + 189*pow(P0x,2.0)*P2x*Q3x*pow(sp,2.0) + 294*pow(P0x,2.0)*P3x*Q2x*pow(sp,2.0) - 189*pow(P1x,2.0)*P2x*Q2x*pow(sp,2.0) + 42*P0x*P1x*pow(Q3x,2.0)*pow(sp,4.0) -
1706 147*P0x*pow(P1x,2.0)*Q3x*pow(sp,4.0) + 42*pow(P0x,2.0)*P1x*Q3x*pow(sp,4.0) - 210*pow(P0x,2.0)*P2x*Q2x*pow(sp,4.0) - 420*pow(P0x,2.0)*P2x*Q3x*pow(sp,3.0) - 210*pow(P0x,2.0)*P3x*Q2x*pow(sp,3.0) -
1707 294*pow(P0x,2.0)*P3x*Q3x*pow(sp,2.0) + 420*pow(P1x,2.0)*P2x*Q2x*pow(sp,3.0) + 189*pow(P1x,2.0)*P2x*Q3x*pow(sp,2.0) + 294*pow(P1x,2.0)*P3x*Q2x*pow(sp,2.0) + 210*pow(P0x,2.0)*P2x*Q3x*pow(sp,4.0) +
1708 210*pow(P0x,2.0)*P3x*Q3x*pow(sp,3.0) - 210*pow(P1x,2.0)*P2x*Q2x*pow(sp,4.0) - 420*pow(P1x,2.0)*P2x*Q3x*pow(sp,3.0) - 210*pow(P1x,2.0)*P3x*Q2x*pow(sp,3.0) - 294*pow(P1x,2.0)*P3x*Q3x*pow(sp,2.0) +
1709 210*pow(P1x,2.0)*P2x*Q3x*pow(sp,4.0) + 210*pow(P1x,2.0)*P3x*Q3x*pow(sp,3.0) - 63*P0x*pow(P0y,2.0)*Q2x*pow(sp,2.0) + 96*P0x*P1x*pow(Q2y,2.0)*pow(sp,2.0) + 70*P0x*pow(P0y,2.0)*Q2x*pow(sp,3.0) +
1710 63*P0x*pow(P0y,2.0)*Q3x*pow(sp,2.0) + 126*P0x*pow(P1y,2.0)*Q2x*pow(sp,2.0) + 63*P1x*pow(P0y,2.0)*Q2x*pow(sp,2.0) + 48*P0x*P1x*pow(Q2y,2.0)*pow(sp,3.0) - 24*P0x*P1x*pow(Q3y,2.0)*pow(sp,2.0) -
1711 21*P0x*pow(P0y,2.0)*Q2x*pow(sp,4.0) - 70*P0x*pow(P0y,2.0)*Q3x*pow(sp,3.0) - 98*P0x*pow(P1y,2.0)*Q2x*pow(sp,3.0) - 126*P0x*pow(P1y,2.0)*Q3x*pow(sp,2.0) + 42*P1x*pow(P0y,2.0)*Q2x*pow(sp,3.0) -
1712 63*P1x*pow(P0y,2.0)*Q3x*pow(sp,2.0) - 126*P1x*pow(P1y,2.0)*Q2x*pow(sp,2.0) - 96*P0x*P1x*pow(Q2y,2.0)*pow(sp,4.0) + 16*P0x*P1x*pow(Q3y,2.0)*pow(sp,3.0) + 21*P0x*pow(P0y,2.0)*Q3x*pow(sp,4.0) +
1713 24*P0x*pow(P1y,2.0)*Q2x*pow(sp,4.0) + 98*P0x*pow(P1y,2.0)*Q3x*pow(sp,3.0) - 39*P1x*pow(P0y,2.0)*Q2x*pow(sp,4.0) - 42*P1x*pow(P0y,2.0)*Q3x*pow(sp,3.0) + 210*P1x*pow(P1y,2.0)*Q2x*pow(sp,3.0) +
1714 126*P1x*pow(P1y,2.0)*Q3x*pow(sp,2.0) + 168*P2x*pow(P0y,2.0)*Q2x*pow(sp,3.0) + 96*P3x*pow(P0y,2.0)*Q2x*pow(sp,2.0) + 24*P0x*P1x*pow(Q3y,2.0)*pow(sp,4.0) - 24*P0x*pow(P1y,2.0)*Q3x*pow(sp,4.0) +
1715 39*P1x*pow(P0y,2.0)*Q3x*pow(sp,4.0) - 84*P1x*pow(P1y,2.0)*Q2x*pow(sp,4.0) - 210*P1x*pow(P1y,2.0)*Q3x*pow(sp,3.0) - 120*P2x*pow(P0y,2.0)*Q2x*pow(sp,4.0) - 168*P2x*pow(P0y,2.0)*Q3x*pow(sp,3.0) +
1716 168*P2x*pow(P1y,2.0)*Q2x*pow(sp,3.0) - 120*P3x*pow(P0y,2.0)*Q2x*pow(sp,3.0) - 96*P3x*pow(P0y,2.0)*Q3x*pow(sp,2.0) + 96*P3x*pow(P1y,2.0)*Q2x*pow(sp,2.0) + 84*P1x*pow(P1y,2.0)*Q3x*pow(sp,4.0) +
1717 120*P2x*pow(P0y,2.0)*Q3x*pow(sp,4.0) - 120*P2x*pow(P1y,2.0)*Q2x*pow(sp,4.0) - 168*P2x*pow(P1y,2.0)*Q3x*pow(sp,3.0) + 120*P3x*pow(P0y,2.0)*Q3x*pow(sp,3.0) - 120*P3x*pow(P1y,2.0)*Q2x*pow(sp,3.0) -
1718 96*P3x*pow(P1y,2.0)*Q3x*pow(sp,2.0) + 120*P2x*pow(P1y,2.0)*Q3x*pow(sp,4.0) + 120*P3x*pow(P1y,2.0)*Q3x*pow(sp,3.0) - 63*pow(P0x,2.0)*P0y*Q2y*pow(sp,2.0) + 96*P0y*P1y*pow(Q2x,2.0)*pow(sp,2.0) +
1719 70*pow(P0x,2.0)*P0y*Q2y*pow(sp,3.0) + 63*pow(P0x,2.0)*P0y*Q3y*pow(sp,2.0) + 63*pow(P0x,2.0)*P1y*Q2y*pow(sp,2.0) + 126*pow(P1x,2.0)*P0y*Q2y*pow(sp,2.0) + 48*P0y*P1y*pow(Q2x,2.0)*pow(sp,3.0) -
1720 24*P0y*P1y*pow(Q3x,2.0)*pow(sp,2.0) - 21*pow(P0x,2.0)*P0y*Q2y*pow(sp,4.0) - 70*pow(P0x,2.0)*P0y*Q3y*pow(sp,3.0) + 42*pow(P0x,2.0)*P1y*Q2y*pow(sp,3.0) - 63*pow(P0x,2.0)*P1y*Q3y*pow(sp,2.0) -
1721 98*pow(P1x,2.0)*P0y*Q2y*pow(sp,3.0) - 126*pow(P1x,2.0)*P0y*Q3y*pow(sp,2.0) - 126*pow(P1x,2.0)*P1y*Q2y*pow(sp,2.0) - 96*P0y*P1y*pow(Q2x,2.0)*pow(sp,4.0) + 16*P0y*P1y*pow(Q3x,2.0)*pow(sp,3.0) +
1722 21*pow(P0x,2.0)*P0y*Q3y*pow(sp,4.0) - 39*pow(P0x,2.0)*P1y*Q2y*pow(sp,4.0) - 42*pow(P0x,2.0)*P1y*Q3y*pow(sp,3.0) + 168*pow(P0x,2.0)*P2y*Q2y*pow(sp,3.0) + 96*pow(P0x,2.0)*P3y*Q2y*pow(sp,2.0) +
1723 24*pow(P1x,2.0)*P0y*Q2y*pow(sp,4.0) + 98*pow(P1x,2.0)*P0y*Q3y*pow(sp,3.0) + 210*pow(P1x,2.0)*P1y*Q2y*pow(sp,3.0) + 126*pow(P1x,2.0)*P1y*Q3y*pow(sp,2.0) + 24*P0y*P1y*pow(Q3x,2.0)*pow(sp,4.0) +
1724 39*pow(P0x,2.0)*P1y*Q3y*pow(sp,4.0) - 120*pow(P0x,2.0)*P2y*Q2y*pow(sp,4.0) - 168*pow(P0x,2.0)*P2y*Q3y*pow(sp,3.0) - 120*pow(P0x,2.0)*P3y*Q2y*pow(sp,3.0) - 96*pow(P0x,2.0)*P3y*Q3y*pow(sp,2.0) -
1725 24*pow(P1x,2.0)*P0y*Q3y*pow(sp,4.0) - 84*pow(P1x,2.0)*P1y*Q2y*pow(sp,4.0) - 210*pow(P1x,2.0)*P1y*Q3y*pow(sp,3.0) + 168*pow(P1x,2.0)*P2y*Q2y*pow(sp,3.0) + 96*pow(P1x,2.0)*P3y*Q2y*pow(sp,2.0) +
1726 120*pow(P0x,2.0)*P2y*Q3y*pow(sp,4.0) + 120*pow(P0x,2.0)*P3y*Q3y*pow(sp,3.0) + 84*pow(P1x,2.0)*P1y*Q3y*pow(sp,4.0) - 120*pow(P1x,2.0)*P2y*Q2y*pow(sp,4.0) - 168*pow(P1x,2.0)*P2y*Q3y*pow(sp,3.0) -
1727 120*pow(P1x,2.0)*P3y*Q2y*pow(sp,3.0) - 96*pow(P1x,2.0)*P3y*Q3y*pow(sp,2.0) + 120*pow(P1x,2.0)*P2y*Q3y*pow(sp,4.0) + 120*pow(P1x,2.0)*P3y*Q3y*pow(sp,3.0) + 42*P0y*P1y*pow(Q2y,2.0)*pow(sp,2.0) +
1728 189*P0y*pow(P1y,2.0)*Q2y*pow(sp,2.0) + 210*P0y*P1y*pow(Q2y,2.0)*pow(sp,3.0) - 42*P0y*P1y*pow(Q3y,2.0)*pow(sp,2.0) - 350*P0y*pow(P1y,2.0)*Q2y*pow(sp,3.0) - 189*P0y*pow(P1y,2.0)*Q3y*pow(sp,2.0) +
1729 70*pow(P0y,2.0)*P1y*Q2y*pow(sp,3.0) - 189*pow(P0y,2.0)*P2y*Q2y*pow(sp,2.0) - 168*P0y*P1y*pow(Q2y,2.0)*pow(sp,4.0) - 14*P0y*P1y*pow(Q3y,2.0)*pow(sp,3.0) + 147*P0y*pow(P1y,2.0)*Q2y*pow(sp,4.0) +
1730 350*P0y*pow(P1y,2.0)*Q3y*pow(sp,3.0) - 42*pow(P0y,2.0)*P1y*Q2y*pow(sp,4.0) - 70*pow(P0y,2.0)*P1y*Q3y*pow(sp,3.0) + 420*pow(P0y,2.0)*P2y*Q2y*pow(sp,3.0) + 189*pow(P0y,2.0)*P2y*Q3y*pow(sp,2.0) +
1731 294*pow(P0y,2.0)*P3y*Q2y*pow(sp,2.0) - 189*pow(P1y,2.0)*P2y*Q2y*pow(sp,2.0) + 42*P0y*P1y*pow(Q3y,2.0)*pow(sp,4.0) - 147*P0y*pow(P1y,2.0)*Q3y*pow(sp,4.0) + 42*pow(P0y,2.0)*P1y*Q3y*pow(sp,4.0) -
1732 210*pow(P0y,2.0)*P2y*Q2y*pow(sp,4.0) - 420*pow(P0y,2.0)*P2y*Q3y*pow(sp,3.0) - 210*pow(P0y,2.0)*P3y*Q2y*pow(sp,3.0) - 294*pow(P0y,2.0)*P3y*Q3y*pow(sp,2.0) + 420*pow(P1y,2.0)*P2y*Q2y*pow(sp,3.0) +
1733 189*pow(P1y,2.0)*P2y*Q3y*pow(sp,2.0) + 294*pow(P1y,2.0)*P3y*Q2y*pow(sp,2.0) + 210*pow(P0y,2.0)*P2y*Q3y*pow(sp,4.0) + 210*pow(P0y,2.0)*P3y*Q3y*pow(sp,3.0) - 210*pow(P1y,2.0)*P2y*Q2y*pow(sp,4.0) -
1734 420*pow(P1y,2.0)*P2y*Q3y*pow(sp,3.0) - 210*pow(P1y,2.0)*P3y*Q2y*pow(sp,3.0) - 294*pow(P1y,2.0)*P3y*Q3y*pow(sp,2.0) + 210*pow(P1y,2.0)*P2y*Q3y*pow(sp,4.0) + 210*pow(P1y,2.0)*P3y*Q3y*pow(sp,3.0) +
1735 126*pow(P0x,2.0)*Q1x*Q2x*pow(sp,2.0) - 378*pow(P0x,2.0)*Q1x*Q2x*pow(sp,3.0) - 126*pow(P0x,2.0)*Q1x*Q3x*pow(sp,2.0) + 126*pow(P1x,2.0)*Q1x*Q2x*pow(sp,2.0) + 210*pow(P0x,2.0)*Q1x*Q2x*pow(sp,4.0) +
1736 378*pow(P0x,2.0)*Q1x*Q3x*pow(sp,3.0) - 378*pow(P1x,2.0)*Q1x*Q2x*pow(sp,3.0) - 126*pow(P1x,2.0)*Q1x*Q3x*pow(sp,2.0) - 210*pow(P0x,2.0)*Q1x*Q3x*pow(sp,4.0) + 98*pow(P0x,2.0)*Q2x*Q3x*pow(sp,3.0) +
1737 210*pow(P1x,2.0)*Q1x*Q2x*pow(sp,4.0) + 378*pow(P1x,2.0)*Q1x*Q3x*pow(sp,3.0) - 63*pow(P0x,2.0)*Q2x*Q3x*pow(sp,4.0) - 210*pow(P1x,2.0)*Q1x*Q3x*pow(sp,4.0) + 98*pow(P1x,2.0)*Q2x*Q3x*pow(sp,3.0) -
1738 63*pow(P1x,2.0)*Q2x*Q3x*pow(sp,4.0) - 36*pow(P0y,2.0)*Q1x*Q2x*pow(sp,2.0) - 144*pow(P0y,2.0)*Q1x*Q2x*pow(sp,3.0) + 36*pow(P0y,2.0)*Q1x*Q3x*pow(sp,2.0) - 36*pow(P1y,2.0)*Q1x*Q2x*pow(sp,2.0) +
1739 120*pow(P0y,2.0)*Q1x*Q2x*pow(sp,4.0) + 144*pow(P0y,2.0)*Q1x*Q3x*pow(sp,3.0) + 36*pow(P0y,2.0)*Q2x*Q3x*pow(sp,2.0) - 144*pow(P1y,2.0)*Q1x*Q2x*pow(sp,3.0) + 36*pow(P1y,2.0)*Q1x*Q3x*pow(sp,2.0) -
1740 120*pow(P0y,2.0)*Q1x*Q3x*pow(sp,4.0) + 32*pow(P0y,2.0)*Q2x*Q3x*pow(sp,3.0) + 120*pow(P1y,2.0)*Q1x*Q2x*pow(sp,4.0) + 144*pow(P1y,2.0)*Q1x*Q3x*pow(sp,3.0) + 36*pow(P1y,2.0)*Q2x*Q3x*pow(sp,2.0) -
1741 36*pow(P0y,2.0)*Q2x*Q3x*pow(sp,4.0) - 120*pow(P1y,2.0)*Q1x*Q3x*pow(sp,4.0) + 32*pow(P1y,2.0)*Q2x*Q3x*pow(sp,3.0) - 36*pow(P1y,2.0)*Q2x*Q3x*pow(sp,4.0) - 36*pow(P0x,2.0)*Q1y*Q2y*pow(sp,2.0) -
1742 144*pow(P0x,2.0)*Q1y*Q2y*pow(sp,3.0) + 36*pow(P0x,2.0)*Q1y*Q3y*pow(sp,2.0) - 36*pow(P1x,2.0)*Q1y*Q2y*pow(sp,2.0) + 120*pow(P0x,2.0)*Q1y*Q2y*pow(sp,4.0) + 144*pow(P0x,2.0)*Q1y*Q3y*pow(sp,3.0) +
1743 36*pow(P0x,2.0)*Q2y*Q3y*pow(sp,2.0) - 144*pow(P1x,2.0)*Q1y*Q2y*pow(sp,3.0) + 36*pow(P1x,2.0)*Q1y*Q3y*pow(sp,2.0) - 120*pow(P0x,2.0)*Q1y*Q3y*pow(sp,4.0) + 32*pow(P0x,2.0)*Q2y*Q3y*pow(sp,3.0) +
1744 120*pow(P1x,2.0)*Q1y*Q2y*pow(sp,4.0) + 144*pow(P1x,2.0)*Q1y*Q3y*pow(sp,3.0) + 36*pow(P1x,2.0)*Q2y*Q3y*pow(sp,2.0) - 36*pow(P0x,2.0)*Q2y*Q3y*pow(sp,4.0) - 120*pow(P1x,2.0)*Q1y*Q3y*pow(sp,4.0) +
1745 32*pow(P1x,2.0)*Q2y*Q3y*pow(sp,3.0) - 36*pow(P1x,2.0)*Q2y*Q3y*pow(sp,4.0) + 126*pow(P0y,2.0)*Q1y*Q2y*pow(sp,2.0) - 378*pow(P0y,2.0)*Q1y*Q2y*pow(sp,3.0) - 126*pow(P0y,2.0)*Q1y*Q3y*pow(sp,2.0) +
1746 126*pow(P1y,2.0)*Q1y*Q2y*pow(sp,2.0) + 210*pow(P0y,2.0)*Q1y*Q2y*pow(sp,4.0) + 378*pow(P0y,2.0)*Q1y*Q3y*pow(sp,3.0) - 378*pow(P1y,2.0)*Q1y*Q2y*pow(sp,3.0) - 126*pow(P1y,2.0)*Q1y*Q3y*pow(sp,2.0) -
1747 210*pow(P0y,2.0)*Q1y*Q3y*pow(sp,4.0) + 98*pow(P0y,2.0)*Q2y*Q3y*pow(sp,3.0) + 210*pow(P1y,2.0)*Q1y*Q2y*pow(sp,4.0) + 378*pow(P1y,2.0)*Q1y*Q3y*pow(sp,3.0) - 63*pow(P0y,2.0)*Q2y*Q3y*pow(sp,4.0) -
1748 210*pow(P1y,2.0)*Q1y*Q3y*pow(sp,4.0) + 98*pow(P1y,2.0)*Q2y*Q3y*pow(sp,3.0) - 63*pow(P1y,2.0)*Q2y*Q3y*pow(sp,4.0) + 28*P0x*P1x*P3x*Q2x - 28*P0x*P1x*P3x*Q3x + 38*P0x*P1x*P0y*Q2y + 38*P0x*P0y*P1y*Q2x -
1749 38*P0x*P1x*P0y*Q3y + 66*P0x*P1x*P1y*Q2y - 38*P0x*P0y*P1y*Q3x + 66*P1x*P0y*P1y*Q2x - 66*P0x*P1x*P1y*Q3y - 30*P0x*P3x*P0y*Q2y - 30*P0x*P0y*P3y*Q2x - 66*P1x*P0y*P1y*Q3x - 32*P0x*P1x*P3y*Q2y + 30*P0x*P3x*P0y*Q3y +
1750 30*P0x*P3x*P1y*Q2y + 30*P0x*P0y*P3y*Q3x + 30*P0x*P1y*P3y*Q2x + 30*P1x*P3x*P0y*Q2y + 30*P1x*P0y*P3y*Q2x - 32*P3x*P0y*P1y*Q2x + 32*P0x*P1x*P3y*Q3y - 30*P0x*P3x*P1y*Q3y - 30*P0x*P1y*P3y*Q3x - 30*P1x*P3x*P0y*Q3y -
1751 30*P1x*P3x*P1y*Q2y - 30*P1x*P0y*P3y*Q3x - 30*P1x*P1y*P3y*Q2x + 32*P3x*P0y*P1y*Q3x + 30*P1x*P3x*P1y*Q3y + 30*P1x*P1y*P3y*Q3x + 28*P0y*P1y*P3y*Q2y - 28*P0y*P1y*P3y*Q3y + 84*P0x*P1x*Q2x*Q3x - 36*P0x*P0y*Q1x*Q2y -
1752 36*P0x*P0y*Q2x*Q1y - 72*P0x*P1x*Q1y*Q2y + 36*P0x*P0y*Q1x*Q3y - 54*P0x*P0y*Q2x*Q2y + 36*P0x*P0y*Q3x*Q1y + 36*P0x*P1y*Q1x*Q2y + 36*P0x*P1y*Q2x*Q1y + 36*P1x*P0y*Q1x*Q2y + 36*P1x*P0y*Q2x*Q1y - 72*P0y*P1y*Q1x*Q2x +
1753 72*P0x*P1x*Q1y*Q3y + 54*P0x*P0y*Q2x*Q3y + 54*P0x*P0y*Q3x*Q2y - 36*P0x*P1y*Q1x*Q3y + 54*P0x*P1y*Q2x*Q2y - 36*P0x*P1y*Q3x*Q1y - 36*P1x*P0y*Q1x*Q3y + 54*P1x*P0y*Q2x*Q2y - 36*P1x*P0y*Q3x*Q1y - 36*P1x*P1y*Q1x*Q2y -
1754 36*P1x*P1y*Q2x*Q1y + 72*P0y*P1y*Q1x*Q3x + 192*P0x*P1x*Q2y*Q3y - 54*P0x*P0y*Q3x*Q3y - 54*P0x*P1y*Q2x*Q3y - 54*P0x*P1y*Q3x*Q2y - 54*P1x*P0y*Q2x*Q3y - 54*P1x*P0y*Q3x*Q2y + 36*P1x*P1y*Q1x*Q3y - 54*P1x*P1y*Q2x*Q2y +
1755 36*P1x*P1y*Q3x*Q1y + 192*P0y*P1y*Q2x*Q3x + 54*P0x*P1y*Q3x*Q3y + 54*P1x*P0y*Q3x*Q3y + 54*P1x*P1y*Q2x*Q3y + 54*P1x*P1y*Q3x*Q2y - 54*P1x*P1y*Q3x*Q3y + 84*P0y*P1y*Q2y*Q3y + 84*P0x*P1x*P3x*Q2x*sp - 84*P0x*P1x*P3x*Q3x*sp -
1756 90*P0x*P3x*P0y*Q2y*sp - 90*P0x*P0y*P3y*Q2x*sp - 96*P0x*P1x*P3y*Q2y*sp + 90*P0x*P3x*P0y*Q3y*sp + 90*P0x*P3x*P1y*Q2y*sp + 90*P0x*P0y*P3y*Q3x*sp + 90*P0x*P1y*P3y*Q2x*sp + 90*P1x*P3x*P0y*Q2y*sp + 90*P1x*P0y*P3y*Q2x*sp -
1757 96*P3x*P0y*P1y*Q2x*sp + 96*P0x*P1x*P3y*Q3y*sp - 90*P0x*P3x*P1y*Q3y*sp - 90*P0x*P1y*P3y*Q3x*sp - 90*P1x*P3x*P0y*Q3y*sp - 90*P1x*P3x*P1y*Q2y*sp - 90*P1x*P0y*P3y*Q3x*sp - 90*P1x*P1y*P3y*Q2x*sp + 96*P3x*P0y*P1y*Q3x*sp +
1758 90*P1x*P3x*P1y*Q3y*sp + 90*P1x*P1y*P3y*Q3x*sp + 84*P0y*P1y*P3y*Q2y*sp - 84*P0y*P1y*P3y*Q3y*sp - 84*P0x*P1x*Q1x*Q2x*sp + 84*P0x*P1x*Q1x*Q3x*sp + 84*P0x*P1x*Q2x*Q3x*sp + 18*P0x*P0y*Q1x*Q2y*sp + 18*P0x*P0y*Q2x*Q1y*sp -
1759 48*P0x*P1x*Q1y*Q2y*sp - 18*P0x*P0y*Q1x*Q3y*sp + 90*P0x*P0y*Q2x*Q2y*sp - 18*P0x*P0y*Q3x*Q1y*sp - 18*P0x*P1y*Q1x*Q2y*sp - 18*P0x*P1y*Q2x*Q1y*sp - 18*P1x*P0y*Q1x*Q2y*sp - 18*P1x*P0y*Q2x*Q1y*sp - 48*P0y*P1y*Q1x*Q2x*sp +
1760 48*P0x*P1x*Q1y*Q3y*sp - 18*P0x*P0y*Q2x*Q3y*sp - 18*P0x*P0y*Q3x*Q2y*sp + 18*P0x*P1y*Q1x*Q3y*sp - 90*P0x*P1y*Q2x*Q2y*sp + 18*P0x*P1y*Q3x*Q1y*sp + 18*P1x*P0y*Q1x*Q3y*sp - 90*P1x*P0y*Q2x*Q2y*sp + 18*P1x*P0y*Q3x*Q1y*sp +
1761 18*P1x*P1y*Q1x*Q2y*sp + 18*P1x*P1y*Q2x*Q1y*sp + 48*P0y*P1y*Q1x*Q3x*sp + 48*P0x*P1x*Q2y*Q3y*sp - 54*P0x*P0y*Q3x*Q3y*sp + 18*P0x*P1y*Q2x*Q3y*sp + 18*P0x*P1y*Q3x*Q2y*sp + 18*P1x*P0y*Q2x*Q3y*sp + 18*P1x*P0y*Q3x*Q2y*sp -
1762 18*P1x*P1y*Q1x*Q3y*sp + 90*P1x*P1y*Q2x*Q2y*sp - 18*P1x*P1y*Q3x*Q1y*sp + 48*P0y*P1y*Q2x*Q3x*sp + 54*P0x*P1y*Q3x*Q3y*sp + 54*P1x*P0y*Q3x*Q3y*sp - 18*P1x*P1y*Q2x*Q3y*sp - 18*P1x*P1y*Q3x*Q2y*sp - 54*P1x*P1y*Q3x*Q3y*sp -
1763 84*P0y*P1y*Q1y*Q2y*sp + 84*P0y*P1y*Q1y*Q3y*sp + 84*P0y*P1y*Q2y*Q3y*sp + 378*P0x*P1x*P2x*Q2x*pow(sp,2.0) - 840*P0x*P1x*P2x*Q2x*pow(sp,3.0) - 378*P0x*P1x*P2x*Q3x*pow(sp,2.0) - 588*P0x*P1x*P3x*Q2x*pow(sp,2.0) +
1764 420*P0x*P1x*P2x*Q2x*pow(sp,4.0) + 840*P0x*P1x*P2x*Q3x*pow(sp,3.0) + 420*P0x*P1x*P3x*Q2x*pow(sp,3.0) + 588*P0x*P1x*P3x*Q3x*pow(sp,2.0) - 420*P0x*P1x*P2x*Q3x*pow(sp,4.0) - 420*P0x*P1x*P3x*Q3x*pow(sp,3.0) -
1765 63*P0x*P1x*P0y*Q2y*pow(sp,2.0) - 63*P0x*P0y*P1y*Q2x*pow(sp,2.0) + 28*P0x*P1x*P0y*Q2y*pow(sp,3.0) + 63*P0x*P1x*P0y*Q3y*pow(sp,2.0) + 63*P0x*P1x*P1y*Q2y*pow(sp,2.0) - 189*P0x*P2x*P0y*Q2y*pow(sp,2.0) +
1766 28*P0x*P0y*P1y*Q2x*pow(sp,3.0) + 63*P0x*P0y*P1y*Q3x*pow(sp,2.0) - 189*P0x*P0y*P2y*Q2x*pow(sp,2.0) + 63*P1x*P0y*P1y*Q2x*pow(sp,2.0) - 3*P0x*P1x*P0y*Q2y*pow(sp,4.0) - 28*P0x*P1x*P0y*Q3y*pow(sp,3.0) -
1767 252*P0x*P1x*P1y*Q2y*pow(sp,3.0) - 63*P0x*P1x*P1y*Q3y*pow(sp,2.0) + 252*P0x*P2x*P0y*Q2y*pow(sp,3.0) + 189*P0x*P2x*P0y*Q3y*pow(sp,2.0) + 189*P0x*P2x*P1y*Q2y*pow(sp,2.0) + 198*P0x*P3x*P0y*Q2y*pow(sp,2.0) -
1768 3*P0x*P0y*P1y*Q2x*pow(sp,4.0) - 28*P0x*P0y*P1y*Q3x*pow(sp,3.0) + 252*P0x*P0y*P2y*Q2x*pow(sp,3.0) + 189*P0x*P0y*P2y*Q3x*pow(sp,2.0) + 198*P0x*P0y*P3y*Q2x*pow(sp,2.0) + 189*P0x*P1y*P2y*Q2x*pow(sp,2.0) +
1769 189*P1x*P2x*P0y*Q2y*pow(sp,2.0) - 252*P1x*P0y*P1y*Q2x*pow(sp,3.0) - 63*P1x*P0y*P1y*Q3x*pow(sp,2.0) + 189*P1x*P0y*P2y*Q2x*pow(sp,2.0) + 3*P0x*P1x*P0y*Q3y*pow(sp,4.0) + 123*P0x*P1x*P1y*Q2y*pow(sp,4.0) +
1770 252*P0x*P1x*P1y*Q3y*pow(sp,3.0) - 336*P0x*P1x*P2y*Q2y*pow(sp,3.0) - 192*P0x*P1x*P3y*Q2y*pow(sp,2.0) - 90*P0x*P2x*P0y*Q2y*pow(sp,4.0) - 252*P0x*P2x*P0y*Q3y*pow(sp,3.0) - 252*P0x*P2x*P1y*Q2y*pow(sp,3.0) -
1771 189*P0x*P2x*P1y*Q3y*pow(sp,2.0) - 90*P0x*P3x*P0y*Q2y*pow(sp,3.0) - 198*P0x*P3x*P0y*Q3y*pow(sp,2.0) - 198*P0x*P3x*P1y*Q2y*pow(sp,2.0) + 3*P0x*P0y*P1y*Q3x*pow(sp,4.0) - 90*P0x*P0y*P2y*Q2x*pow(sp,4.0) -
1772 252*P0x*P0y*P2y*Q3x*pow(sp,3.0) - 90*P0x*P0y*P3y*Q2x*pow(sp,3.0) - 198*P0x*P0y*P3y*Q3x*pow(sp,2.0) - 252*P0x*P1y*P2y*Q2x*pow(sp,3.0) - 189*P0x*P1y*P2y*Q3x*pow(sp,2.0) - 198*P0x*P1y*P3y*Q2x*pow(sp,2.0) -
1773 252*P1x*P2x*P0y*Q2y*pow(sp,3.0) - 189*P1x*P2x*P0y*Q3y*pow(sp,2.0) - 189*P1x*P2x*P1y*Q2y*pow(sp,2.0) - 198*P1x*P3x*P0y*Q2y*pow(sp,2.0) + 123*P1x*P0y*P1y*Q2x*pow(sp,4.0) + 252*P1x*P0y*P1y*Q3x*pow(sp,3.0) -
1774 252*P1x*P0y*P2y*Q2x*pow(sp,3.0) - 189*P1x*P0y*P2y*Q3x*pow(sp,2.0) - 198*P1x*P0y*P3y*Q2x*pow(sp,2.0) - 189*P1x*P1y*P2y*Q2x*pow(sp,2.0) - 336*P2x*P0y*P1y*Q2x*pow(sp,3.0) - 192*P3x*P0y*P1y*Q2x*pow(sp,2.0) -
1775 123*P0x*P1x*P1y*Q3y*pow(sp,4.0) + 240*P0x*P1x*P2y*Q2y*pow(sp,4.0) + 336*P0x*P1x*P2y*Q3y*pow(sp,3.0) + 240*P0x*P1x*P3y*Q2y*pow(sp,3.0) + 192*P0x*P1x*P3y*Q3y*pow(sp,2.0) + 90*P0x*P2x*P0y*Q3y*pow(sp,4.0) +
1776 90*P0x*P2x*P1y*Q2y*pow(sp,4.0) + 252*P0x*P2x*P1y*Q3y*pow(sp,3.0) + 90*P0x*P3x*P0y*Q3y*pow(sp,3.0) + 90*P0x*P3x*P1y*Q2y*pow(sp,3.0) + 198*P0x*P3x*P1y*Q3y*pow(sp,2.0) + 90*P0x*P0y*P2y*Q3x*pow(sp,4.0) +
1777 90*P0x*P0y*P3y*Q3x*pow(sp,3.0) + 90*P0x*P1y*P2y*Q2x*pow(sp,4.0) + 252*P0x*P1y*P2y*Q3x*pow(sp,3.0) + 90*P0x*P1y*P3y*Q2x*pow(sp,3.0) + 198*P0x*P1y*P3y*Q3x*pow(sp,2.0) + 90*P1x*P2x*P0y*Q2y*pow(sp,4.0) +
1778 252*P1x*P2x*P0y*Q3y*pow(sp,3.0) + 252*P1x*P2x*P1y*Q2y*pow(sp,3.0) + 189*P1x*P2x*P1y*Q3y*pow(sp,2.0) + 90*P1x*P3x*P0y*Q2y*pow(sp,3.0) + 198*P1x*P3x*P0y*Q3y*pow(sp,2.0) + 198*P1x*P3x*P1y*Q2y*pow(sp,2.0) -
1779 123*P1x*P0y*P1y*Q3x*pow(sp,4.0) + 90*P1x*P0y*P2y*Q2x*pow(sp,4.0) + 252*P1x*P0y*P2y*Q3x*pow(sp,3.0) + 90*P1x*P0y*P3y*Q2x*pow(sp,3.0) + 198*P1x*P0y*P3y*Q3x*pow(sp,2.0) + 252*P1x*P1y*P2y*Q2x*pow(sp,3.0) +
1780 189*P1x*P1y*P2y*Q3x*pow(sp,2.0) + 198*P1x*P1y*P3y*Q2x*pow(sp,2.0) + 240*P2x*P0y*P1y*Q2x*pow(sp,4.0) + 336*P2x*P0y*P1y*Q3x*pow(sp,3.0) + 240*P3x*P0y*P1y*Q2x*pow(sp,3.0) + 192*P3x*P0y*P1y*Q3x*pow(sp,2.0) -
1781 240*P0x*P1x*P2y*Q3y*pow(sp,4.0) - 240*P0x*P1x*P3y*Q3y*pow(sp,3.0) - 90*P0x*P2x*P1y*Q3y*pow(sp,4.0) - 90*P0x*P3x*P1y*Q3y*pow(sp,3.0) - 90*P0x*P1y*P2y*Q3x*pow(sp,4.0) - 90*P0x*P1y*P3y*Q3x*pow(sp,3.0) -
1782 90*P1x*P2x*P0y*Q3y*pow(sp,4.0) - 90*P1x*P2x*P1y*Q2y*pow(sp,4.0) - 252*P1x*P2x*P1y*Q3y*pow(sp,3.0) - 90*P1x*P3x*P0y*Q3y*pow(sp,3.0) - 90*P1x*P3x*P1y*Q2y*pow(sp,3.0) - 198*P1x*P3x*P1y*Q3y*pow(sp,2.0) -
1783 90*P1x*P0y*P2y*Q3x*pow(sp,4.0) - 90*P1x*P0y*P3y*Q3x*pow(sp,3.0) - 90*P1x*P1y*P2y*Q2x*pow(sp,4.0) - 252*P1x*P1y*P2y*Q3x*pow(sp,3.0) - 90*P1x*P1y*P3y*Q2x*pow(sp,3.0) - 198*P1x*P1y*P3y*Q3x*pow(sp,2.0) -
1784 240*P2x*P0y*P1y*Q3x*pow(sp,4.0) - 240*P3x*P0y*P1y*Q3x*pow(sp,3.0) + 90*P1x*P2x*P1y*Q3y*pow(sp,4.0) + 90*P1x*P3x*P1y*Q3y*pow(sp,3.0) + 90*P1x*P1y*P2y*Q3x*pow(sp,4.0) + 90*P1x*P1y*P3y*Q3x*pow(sp,3.0) +
1785 378*P0y*P1y*P2y*Q2y*pow(sp,2.0) - 840*P0y*P1y*P2y*Q2y*pow(sp,3.0) - 378*P0y*P1y*P2y*Q3y*pow(sp,2.0) - 588*P0y*P1y*P3y*Q2y*pow(sp,2.0) + 420*P0y*P1y*P2y*Q2y*pow(sp,4.0) + 840*P0y*P1y*P2y*Q3y*pow(sp,3.0) +
1786 420*P0y*P1y*P3y*Q2y*pow(sp,3.0) + 588*P0y*P1y*P3y*Q3y*pow(sp,2.0) - 420*P0y*P1y*P2y*Q3y*pow(sp,4.0) - 420*P0y*P1y*P3y*Q3y*pow(sp,3.0) - 252*P0x*P1x*Q1x*Q2x*pow(sp,2.0) + 756*P0x*P1x*Q1x*Q2x*pow(sp,3.0) +
1787 252*P0x*P1x*Q1x*Q3x*pow(sp,2.0) - 420*P0x*P1x*Q1x*Q2x*pow(sp,4.0) - 756*P0x*P1x*Q1x*Q3x*pow(sp,3.0) + 420*P0x*P1x*Q1x*Q3x*pow(sp,4.0) - 196*P0x*P1x*Q2x*Q3x*pow(sp,3.0) + 126*P0x*P1x*Q2x*Q3x*pow(sp,4.0) +
1788 162*P0x*P0y*Q1x*Q2y*pow(sp,2.0) + 162*P0x*P0y*Q2x*Q1y*pow(sp,2.0) + 72*P0x*P1x*Q1y*Q2y*pow(sp,2.0) - 234*P0x*P0y*Q1x*Q2y*pow(sp,3.0) - 162*P0x*P0y*Q1x*Q3y*pow(sp,2.0) - 234*P0x*P0y*Q2x*Q1y*pow(sp,3.0) +
1789 54*P0x*P0y*Q2x*Q2y*pow(sp,2.0) - 162*P0x*P0y*Q3x*Q1y*pow(sp,2.0) - 162*P0x*P1y*Q1x*Q2y*pow(sp,2.0) - 162*P0x*P1y*Q2x*Q1y*pow(sp,2.0) - 162*P1x*P0y*Q1x*Q2y*pow(sp,2.0) - 162*P1x*P0y*Q2x*Q1y*pow(sp,2.0) +
1790 72*P0y*P1y*Q1x*Q2x*pow(sp,2.0) + 288*P0x*P1x*Q1y*Q2y*pow(sp,3.0) - 72*P0x*P1x*Q1y*Q3y*pow(sp,2.0) + 90*P0x*P0y*Q1x*Q2y*pow(sp,4.0) + 234*P0x*P0y*Q1x*Q3y*pow(sp,3.0) + 90*P0x*P0y*Q2x*Q1y*pow(sp,4.0) -
1791 162*P0x*P0y*Q2x*Q2y*pow(sp,3.0) - 36*P0x*P0y*Q2x*Q3y*pow(sp,2.0) + 234*P0x*P0y*Q3x*Q1y*pow(sp,3.0) - 36*P0x*P0y*Q3x*Q2y*pow(sp,2.0) + 234*P0x*P1y*Q1x*Q2y*pow(sp,3.0) + 162*P0x*P1y*Q1x*Q3y*pow(sp,2.0) +
1792 234*P0x*P1y*Q2x*Q1y*pow(sp,3.0) - 54*P0x*P1y*Q2x*Q2y*pow(sp,2.0) + 162*P0x*P1y*Q3x*Q1y*pow(sp,2.0) + 234*P1x*P0y*Q1x*Q2y*pow(sp,3.0) + 162*P1x*P0y*Q1x*Q3y*pow(sp,2.0) + 234*P1x*P0y*Q2x*Q1y*pow(sp,3.0) -
1793 54*P1x*P0y*Q2x*Q2y*pow(sp,2.0) + 162*P1x*P0y*Q3x*Q1y*pow(sp,2.0) + 162*P1x*P1y*Q1x*Q2y*pow(sp,2.0) + 162*P1x*P1y*Q2x*Q1y*pow(sp,2.0) + 288*P0y*P1y*Q1x*Q2x*pow(sp,3.0) - 72*P0y*P1y*Q1x*Q3x*pow(sp,2.0) -
1794 240*P0x*P1x*Q1y*Q2y*pow(sp,4.0) - 288*P0x*P1x*Q1y*Q3y*pow(sp,3.0) - 72*P0x*P1x*Q2y*Q3y*pow(sp,2.0) - 90*P0x*P0y*Q1x*Q3y*pow(sp,4.0) + 72*P0x*P0y*Q2x*Q2y*pow(sp,4.0) + 66*P0x*P0y*Q2x*Q3y*pow(sp,3.0) -
1795 90*P0x*P0y*Q3x*Q1y*pow(sp,4.0) + 66*P0x*P0y*Q3x*Q2y*pow(sp,3.0) + 18*P0x*P0y*Q3x*Q3y*pow(sp,2.0) - 90*P0x*P1y*Q1x*Q2y*pow(sp,4.0) - 234*P0x*P1y*Q1x*Q3y*pow(sp,3.0) - 90*P0x*P1y*Q2x*Q1y*pow(sp,4.0) +
1796 162*P0x*P1y*Q2x*Q2y*pow(sp,3.0) + 36*P0x*P1y*Q2x*Q3y*pow(sp,2.0) - 234*P0x*P1y*Q3x*Q1y*pow(sp,3.0) + 36*P0x*P1y*Q3x*Q2y*pow(sp,2.0) - 90*P1x*P0y*Q1x*Q2y*pow(sp,4.0) - 234*P1x*P0y*Q1x*Q3y*pow(sp,3.0) -
1797 90*P1x*P0y*Q2x*Q1y*pow(sp,4.0) + 162*P1x*P0y*Q2x*Q2y*pow(sp,3.0) + 36*P1x*P0y*Q2x*Q3y*pow(sp,2.0) - 234*P1x*P0y*Q3x*Q1y*pow(sp,3.0) + 36*P1x*P0y*Q3x*Q2y*pow(sp,2.0) - 234*P1x*P1y*Q1x*Q2y*pow(sp,3.0) -
1798 162*P1x*P1y*Q1x*Q3y*pow(sp,2.0) - 234*P1x*P1y*Q2x*Q1y*pow(sp,3.0) + 54*P1x*P1y*Q2x*Q2y*pow(sp,2.0) - 162*P1x*P1y*Q3x*Q1y*pow(sp,2.0) - 240*P0y*P1y*Q1x*Q2x*pow(sp,4.0) - 288*P0y*P1y*Q1x*Q3x*pow(sp,3.0) -
1799 72*P0y*P1y*Q2x*Q3x*pow(sp,2.0) + 240*P0x*P1x*Q1y*Q3y*pow(sp,4.0) - 64*P0x*P1x*Q2y*Q3y*pow(sp,3.0) - 27*P0x*P0y*Q2x*Q3y*pow(sp,4.0) - 27*P0x*P0y*Q3x*Q2y*pow(sp,4.0) + 30*P0x*P0y*Q3x*Q3y*pow(sp,3.0) +
1800 90*P0x*P1y*Q1x*Q3y*pow(sp,4.0) - 72*P0x*P1y*Q2x*Q2y*pow(sp,4.0) - 66*P0x*P1y*Q2x*Q3y*pow(sp,3.0) + 90*P0x*P1y*Q3x*Q1y*pow(sp,4.0) - 66*P0x*P1y*Q3x*Q2y*pow(sp,3.0) - 18*P0x*P1y*Q3x*Q3y*pow(sp,2.0) +
1801 90*P1x*P0y*Q1x*Q3y*pow(sp,4.0) - 72*P1x*P0y*Q2x*Q2y*pow(sp,4.0) - 66*P1x*P0y*Q2x*Q3y*pow(sp,3.0) + 90*P1x*P0y*Q3x*Q1y*pow(sp,4.0) - 66*P1x*P0y*Q3x*Q2y*pow(sp,3.0) - 18*P1x*P0y*Q3x*Q3y*pow(sp,2.0) +
1802 90*P1x*P1y*Q1x*Q2y*pow(sp,4.0) + 234*P1x*P1y*Q1x*Q3y*pow(sp,3.0) + 90*P1x*P1y*Q2x*Q1y*pow(sp,4.0) - 162*P1x*P1y*Q2x*Q2y*pow(sp,3.0) - 36*P1x*P1y*Q2x*Q3y*pow(sp,2.0) + 234*P1x*P1y*Q3x*Q1y*pow(sp,3.0) -
1803 36*P1x*P1y*Q3x*Q2y*pow(sp,2.0) + 240*P0y*P1y*Q1x*Q3x*pow(sp,4.0) - 64*P0y*P1y*Q2x*Q3x*pow(sp,3.0) + 72*P0x*P1x*Q2y*Q3y*pow(sp,4.0) - 18*P0x*P0y*Q3x*Q3y*pow(sp,4.0) + 27*P0x*P1y*Q2x*Q3y*pow(sp,4.0) +
1804 27*P0x*P1y*Q3x*Q2y*pow(sp,4.0) - 30*P0x*P1y*Q3x*Q3y*pow(sp,3.0) + 27*P1x*P0y*Q2x*Q3y*pow(sp,4.0) + 27*P1x*P0y*Q3x*Q2y*pow(sp,4.0) - 30*P1x*P0y*Q3x*Q3y*pow(sp,3.0) - 90*P1x*P1y*Q1x*Q3y*pow(sp,4.0) +
1805 72*P1x*P1y*Q2x*Q2y*pow(sp,4.0) + 66*P1x*P1y*Q2x*Q3y*pow(sp,3.0) - 90*P1x*P1y*Q3x*Q1y*pow(sp,4.0) + 66*P1x*P1y*Q3x*Q2y*pow(sp,3.0) + 18*P1x*P1y*Q3x*Q3y*pow(sp,2.0) + 72*P0y*P1y*Q2x*Q3x*pow(sp,4.0) +
1806 18*P0x*P1y*Q3x*Q3y*pow(sp,4.0) + 18*P1x*P0y*Q3x*Q3y*pow(sp,4.0) - 27*P1x*P1y*Q2x*Q3y*pow(sp,4.0) - 27*P1x*P1y*Q3x*Q2y*pow(sp,4.0) + 30*P1x*P1y*Q3x*Q3y*pow(sp,3.0) - 18*P1x*P1y*Q3x*Q3y*pow(sp,4.0) -
1807 252*P0y*P1y*Q1y*Q2y*pow(sp,2.0) + 756*P0y*P1y*Q1y*Q2y*pow(sp,3.0) + 252*P0y*P1y*Q1y*Q3y*pow(sp,2.0) - 420*P0y*P1y*Q1y*Q2y*pow(sp,4.0) - 756*P0y*P1y*Q1y*Q3y*pow(sp,3.0) + 420*P0y*P1y*Q1y*Q3y*pow(sp,4.0) -
1808 196*P0y*P1y*Q2y*Q3y*pow(sp,3.0) + 126*P0y*P1y*Q2y*Q3y*pow(sp,4.0))/(3*(7*pow(P0x,2.0)*pow(Q2x,2.0) + 7*pow(P0x,2.0)*pow(Q3x,2.0) + 7*pow(P1x,2.0)*pow(Q2x,2.0) + 7*pow(P1x,2.0)*pow(Q3x,2.0) +
1809 16*pow(P0x,2.0)*pow(Q2y,2.0) + 16*pow(P0y,2.0)*pow(Q2x,2.0) + 16*pow(P0x,2.0)*pow(Q3y,2.0) + 16*pow(P1x,2.0)*pow(Q2y,2.0) + 16*pow(P0y,2.0)*pow(Q3x,2.0) + 16*pow(P1y,2.0)*pow(Q2x,2.0) + 16*pow(P1x,2.0)*pow(Q3y,2.0) +
1810 16*pow(P1y,2.0)*pow(Q3x,2.0) + 7*pow(P0y,2.0)*pow(Q2y,2.0) + 7*pow(P0y,2.0)*pow(Q3y,2.0) + 7*pow(P1y,2.0)*pow(Q2y,2.0) + 7*pow(P1y,2.0)*pow(Q3y,2.0) - 14*P0x*P1x*pow(Q2x,2.0) - 14*P0x*P1x*pow(Q3x,2.0) -
1811 32*P0x*P1x*pow(Q2y,2.0) - 32*P0x*P1x*pow(Q3y,2.0) - 32*P0y*P1y*pow(Q2x,2.0) - 32*P0y*P1y*pow(Q3x,2.0) - 14*P0y*P1y*pow(Q2y,2.0) - 14*P0y*P1y*pow(Q3y,2.0) - 14*pow(P0x,2.0)*Q2x*Q3x - 14*pow(P1x,2.0)*Q2x*Q3x -
1812 32*pow(P0y,2.0)*Q2x*Q3x - 32*pow(P1y,2.0)*Q2x*Q3x - 32*pow(P0x,2.0)*Q2y*Q3y - 32*pow(P1x,2.0)*Q2y*Q3y - 14*pow(P0y,2.0)*Q2y*Q3y - 14*pow(P1y,2.0)*Q2y*Q3y + 28*P0x*P1x*Q2x*Q3x - 18*P0x*P0y*Q2x*Q2y + 18*P0x*P0y*Q2x*Q3y +
1813 18*P0x*P0y*Q3x*Q2y + 18*P0x*P1y*Q2x*Q2y + 18*P1x*P0y*Q2x*Q2y + 64*P0x*P1x*Q2y*Q3y - 18*P0x*P0y*Q3x*Q3y - 18*P0x*P1y*Q2x*Q3y - 18*P0x*P1y*Q3x*Q2y - 18*P1x*P0y*Q2x*Q3y - 18*P1x*P0y*Q3x*Q2y - 18*P1x*P1y*Q2x*Q2y +
1814 64*P0y*P1y*Q2x*Q3x + 18*P0x*P1y*Q3x*Q3y + 18*P1x*P0y*Q3x*Q3y + 18*P1x*P1y*Q2x*Q3y + 18*P1x*P1y*Q3x*Q2y - 18*P1x*P1y*Q3x*Q3y + 28*P0y*P1y*Q2y*Q3y));
1815
1816 // Just an arbitrary correction heuristic for negative factors ...
1817 if (out_pfactor < 0)
1818 out_pfactor = -0.1 * out_pfactor;
1819 if (out_qfactor < 0)
1820 out_qfactor = -0.1 * out_qfactor;
1821 }
1822
calcBezierPointDeletionRetainingShapeCost(MapCoord p0,MapCoordF p1,MapCoordF p2,MapCoord p3,PathObject * reference)1823 float PathObject::calcBezierPointDeletionRetainingShapeCost(MapCoord p0, MapCoordF p1, MapCoordF p2, MapCoord p3, PathObject* reference)
1824 {
1825 const int num_test_points = 20;
1826 QBezier curve = QBezier::fromPoints(QPointF(p0), QPointF(p1), QPointF(p2), QPointF(p3));
1827
1828 float cost = 0;
1829 for (int i = 0; i < num_test_points; ++i)
1830 {
1831 auto point = MapCoordF { curve.pointAt((i + 1) / (float)(num_test_points + 1)) };
1832 float distance_sq;
1833 PathCoord path_coord;
1834 reference->calcClosestPointOnPath(MapCoordF(point), distance_sq, path_coord);
1835 cost += distance_sq;
1836 }
1837 // Just some random scaling to pretend that we have 50 sample points
1838 return cost * (50 / 20.0f);
1839 }
1840
calcBezierPointDeletionRetainingShapeOptimization(MapCoord p0,MapCoord p1,MapCoord p2,MapCoord q0,MapCoord q1,MapCoord q2,MapCoord q3,double & out_pfactor,double & out_qfactor)1841 void PathObject::calcBezierPointDeletionRetainingShapeOptimization(MapCoord p0, MapCoord p1, MapCoord p2, MapCoord q0, MapCoord q1, MapCoord q2, MapCoord q3, double& out_pfactor, double& out_qfactor)
1842 {
1843 const float gradient_abort_threshold = 0.05f; // if the gradient magnitude is lower than this over num_abort_steps step, the optimization is aborted
1844 const float decrease_abort_threshold = 0.004f; // if the cost descrease if lower than this over num_abort_steps step, the optimization is aborted
1845 const int num_abort_steps = 2;
1846 const double derivative_delta = 0.05;
1847 const int num_tested_step_sizes = 5;
1848 const int max_num_line_search_iterations = 5;
1849 float step_size_stepping_base = 0.001f;
1850
1851 static LineSymbol line_symbol;
1852 PathObject old_curve(&line_symbol);
1853 p0.setCurveStart(true);
1854 old_curve.addCoordinate(p0);
1855 old_curve.addCoordinate(p1);
1856 old_curve.addCoordinate(p2);
1857 q0.setCurveStart(true);
1858 old_curve.addCoordinate(q0);
1859 old_curve.addCoordinate(q1);
1860 old_curve.addCoordinate(q2);
1861 q3.setCurveStart(false);
1862 q3.setClosePoint(false);
1863 q3.setHolePoint(false);
1864 old_curve.addCoordinate(q3);
1865 old_curve.update();
1866
1867 float cur_cost = 0;
1868 int num_no_improvement_iterations = 0;
1869 float old_gradient[2];
1870 float old_step_dir[2];
1871 for (int i = 0; i < 30; ++i)
1872 {
1873 MapCoordF p_direction = MapCoordF(p1) - MapCoordF(p0);
1874 MapCoordF r1 = MapCoordF(p0) + out_pfactor * p_direction;
1875 MapCoordF q_direction = MapCoordF(q2) - MapCoordF(q3);
1876 MapCoordF r2 = MapCoordF(q3) + out_qfactor * q_direction;
1877
1878 // Calculate gradient and cost (if first iteration)
1879 if (i == 0)
1880 cur_cost = calcBezierPointDeletionRetainingShapeCost(p0, r1, r2, q3, &old_curve);
1881 //if (i == 0)
1882 // qDebug() << "\nStart cost: " << cur_cost;
1883
1884 float gradient[2];
1885 gradient[0] = (calcBezierPointDeletionRetainingShapeCost(p0, r1 + derivative_delta * p_direction, r2, q3, &old_curve) -
1886 calcBezierPointDeletionRetainingShapeCost(p0, r1 - derivative_delta * p_direction, r2, q3, &old_curve)) / (2 * derivative_delta);
1887 gradient[1] = (calcBezierPointDeletionRetainingShapeCost(p0, r1, r2 + derivative_delta * q_direction, q3, &old_curve) -
1888 calcBezierPointDeletionRetainingShapeCost(p0, r1, r2 - derivative_delta * q_direction, q3, &old_curve)) / (2 * derivative_delta);
1889
1890 // Calculate step direction
1891 float step_dir_p;
1892 float step_dir_q;
1893 float conjugate_gradient_factor = 0;
1894 if (i == 0)
1895 {
1896 // Steepest descent
1897 step_dir_p = -gradient[0];
1898 step_dir_q = -gradient[1];
1899 }
1900 else
1901 {
1902 // Conjugate gradient
1903 // Fletcher - Reeves:
1904 //conjugate_gradient_factor = (pow(gradient[0], 2.0) + pow(gradient[1], 2.0)) / (pow(old_gradient[0], 2.0) + pow(old_gradient[1], 2.0));
1905 // Polak – Ribiere:
1906 conjugate_gradient_factor = qMax(0.0, ( gradient[0] * (gradient[0] - old_gradient[0]) + gradient[1] * (gradient[1] - old_gradient[1]) ) /
1907 ( pow(old_gradient[0], 2.0) + pow(old_gradient[1], 2.0) ));
1908 //qDebug() << "Factor: " << conjugate_gradient_factor;
1909 step_dir_p = -gradient[0] + conjugate_gradient_factor * old_step_dir[0];
1910 step_dir_q = -gradient[1] + conjugate_gradient_factor * old_step_dir[1];
1911 }
1912
1913 // Line search in step direction for lowest cost
1914 float best_step_size = 0;
1915 float best_step_size_cost = cur_cost;
1916 int best_step_factor = 0;
1917 float adjusted_step_size_stepping_base = step_size_stepping_base;
1918
1919 for (int iteration = 0; iteration < max_num_line_search_iterations; ++iteration)
1920 {
1921 const float step_size_stepping = adjusted_step_size_stepping_base * (out_pfactor + out_qfactor) / 2; // * qSqrt(step_dir_p*step_dir_p + step_dir_q*step_dir_q); // qSqrt(cur_cost);
1922 for (int step_test = 1; step_test <= num_tested_step_sizes; ++step_test)
1923 {
1924 float step_size = step_test * step_size_stepping;
1925 float test_p_step = step_size * step_dir_p;
1926 float test_q_step = step_size * step_dir_q;
1927
1928 MapCoordF test_r1 = r1 + test_p_step * p_direction;
1929 MapCoordF test_r2 = r2 + test_q_step * q_direction;
1930 float test_cost = calcBezierPointDeletionRetainingShapeCost(p0, test_r1, test_r2, q3, &old_curve);
1931 if (test_cost < best_step_size_cost)
1932 {
1933 best_step_size_cost = test_cost;
1934 best_step_size = step_size;
1935 best_step_factor = step_test;
1936 }
1937 }
1938 //qDebug() << best_step_factor;
1939 if (best_step_factor == num_tested_step_sizes)
1940 adjusted_step_size_stepping_base *= num_tested_step_sizes;
1941 else if (best_step_factor == 0)
1942 adjusted_step_size_stepping_base *= (1 / (float)num_tested_step_sizes);
1943 else
1944 break;
1945 if (iteration < 3)
1946 step_size_stepping_base = adjusted_step_size_stepping_base;
1947 }
1948 if (best_step_factor == 0 && conjugate_gradient_factor == 0)
1949 return;
1950
1951 // Update optimized parameters and constrain them to non-negative values
1952 out_pfactor += best_step_size * step_dir_p;
1953 if (out_pfactor < 0)
1954 out_pfactor = 0;
1955 out_qfactor += best_step_size * step_dir_q;
1956 if (out_qfactor < 0)
1957 out_qfactor = 0;
1958
1959 // Abort if gradient is really low for a number of steps
1960 float gradient_magnitude = qSqrt(gradient[0]*gradient[0] + gradient[1]*gradient[1]);
1961 //qDebug() << "Gradient magnitude: " << gradient_magnitude;
1962 if (gradient_magnitude < gradient_abort_threshold || cur_cost - best_step_size_cost < decrease_abort_threshold)
1963 {
1964 ++num_no_improvement_iterations;
1965 if (num_no_improvement_iterations == num_abort_steps)
1966 break;
1967 }
1968 else
1969 num_no_improvement_iterations = 0;
1970
1971 cur_cost = best_step_size_cost;
1972 old_gradient[0] = gradient[0];
1973 old_gradient[1] = gradient[1];
1974 old_step_dir[0] = step_dir_p;
1975 old_step_dir[1] = step_dir_q;
1976 //qDebug() << "Cost: " << cur_cost;
1977 }
1978 }
1979
appendPath(const PathObject * other)1980 void PathObject::appendPath(const PathObject* other)
1981 {
1982 coords.reserve(coords.size() + other->coords.size());
1983 coords.insert(coords.end(), other->coords.begin(), other->coords.end());
1984
1985 recalculateParts();
1986 setOutputDirty();
1987 }
1988
appendPathPart(const PathPart & part)1989 void PathObject::appendPathPart(const PathPart &part)
1990 {
1991 coords.reserve(coords.size() + part.size());
1992 for (std::size_t i = 0; i < part.size(); ++i)
1993 coords.emplace_back(part.path->coords[part.first_index + i]);
1994
1995 recalculateParts();
1996 setOutputDirty();
1997 }
1998
reverse()1999 void PathObject::reverse()
2000 {
2001 for (auto& part : path_parts)
2002 part.reverse();
2003
2004 Q_ASSERT(isOutputDirty());
2005 }
2006
closeAllParts()2007 void PathObject::closeAllParts()
2008 {
2009 for (auto& part : path_parts)
2010 part.setClosed(true, true);
2011 }
2012
convertToCurves(PathObject ** undo_duplicate)2013 bool PathObject::convertToCurves(PathObject** undo_duplicate)
2014 {
2015 bool converted_a_range = false;
2016 for (const auto& part : path_parts)
2017 {
2018 for (auto index = part.first_index; index < part.last_index; index += 3)
2019 {
2020 if (!coords[index].isCurveStart())
2021 {
2022 auto end_index = index + 1;
2023 while (end_index < part.last_index && !coords[end_index].isCurveStart())
2024 ++end_index;
2025
2026 if (!converted_a_range)
2027 {
2028 if (undo_duplicate)
2029 *undo_duplicate = duplicate()->asPath();
2030 converted_a_range = true;
2031 }
2032
2033 index = convertRangeToCurves(part, index, end_index);
2034 }
2035 }
2036 }
2037
2038 Q_ASSERT(!converted_a_range || isOutputDirty());
2039 return converted_a_range;
2040 }
2041
convertRangeToCurves(const PathPart & part,PathPart::size_type start_index,PathPart::size_type end_index)2042 PathPart::size_type PathObject::convertRangeToCurves(const PathPart& part, PathPart::size_type start_index, PathPart::size_type end_index)
2043 {
2044 Q_ASSERT(end_index > start_index);
2045
2046 Q_ASSERT(!coords[start_index].isCurveStart());
2047 coords[start_index].setCurveStart(true);
2048
2049 // Special case: last coordinate
2050 MapCoord direction;
2051 if (end_index != part.last_index)
2052 {
2053 // Use direction of next segment
2054 direction = coords[end_index + 1] - coords[end_index];
2055 }
2056 else if (part.isClosed())
2057 {
2058 // Use average direction at close point
2059 direction = coords[part.first_index + 1] - coords[end_index - 1];
2060 }
2061 else
2062 {
2063 // Use direction of last segment
2064 direction = coords[end_index] - coords[end_index - 1];
2065 }
2066
2067 MapCoordF tangent = MapCoordF(direction);
2068 tangent.normalize();
2069
2070 MapCoord end_handle = coords[end_index];
2071 end_handle.setFlags(0);
2072 auto baseline = (coords[end_index] - coords[end_index - 1]).length() * BEZIER_HANDLE_DISTANCE;
2073 end_handle = end_handle - MapCoord(tangent * baseline);
2074
2075 // Special case: first coordinate
2076 if (start_index != part.first_index)
2077 {
2078 // Use direction of previous segment
2079 direction = coords[start_index] - coords[start_index - 1];
2080 }
2081 else if (part.isClosed())
2082 {
2083 // Use average direction at close point
2084 direction = coords[start_index + 1] - coords[part.last_index - 1];
2085 }
2086 else
2087 {
2088 // Use direction of first segment
2089 direction = coords[start_index + 1] - coords[start_index];
2090 }
2091
2092 tangent = MapCoordF(direction);
2093 tangent.normalize();
2094
2095 MapCoord handle = coords[start_index];
2096 handle.setFlags(0);
2097 baseline = (coords[start_index + 1] - coords[start_index]).length() * BEZIER_HANDLE_DISTANCE;
2098 handle = handle + MapCoord(tangent * baseline);
2099 addCoordinate(start_index + 1, handle);
2100 ++end_index;
2101
2102 // In-between coordinates
2103 for (MapCoordVector::size_type c = start_index + 2; c < end_index; ++c)
2104 {
2105 direction = coords[c + 1] - coords[c - 2];
2106 tangent = MapCoordF(direction);
2107 tangent.normalize();
2108
2109 // Add previous handle
2110 handle = coords[c];
2111 handle.setFlags(0);
2112 baseline = (coords[c] - coords[c - 2]).length() * BEZIER_HANDLE_DISTANCE;
2113 handle = handle - MapCoord(tangent * baseline);
2114
2115 addCoordinate(c, handle);
2116 ++c;
2117 ++end_index;
2118
2119 Q_ASSERT(!coords[c].isCurveStart());
2120 // Set curve start flag on point
2121 coords[c].setCurveStart(true);
2122
2123 // Add next handle
2124 handle = coords[c];
2125 handle.setFlags(0);
2126 baseline = (coords[c + 1] - coords[c]).length() * BEZIER_HANDLE_DISTANCE;
2127 handle = handle + MapCoord(tangent * baseline);
2128
2129 addCoordinate(c + 1, handle);
2130 ++c;
2131 ++end_index;
2132 }
2133
2134 // Add last handle
2135 addCoordinate(end_index, end_handle);
2136 ++end_index;
2137
2138 return end_index;
2139 }
2140
simplify(PathObject ** undo_duplicate,double threshold)2141 bool PathObject::simplify(PathObject** undo_duplicate, double threshold)
2142 {
2143 // A copy for reference and undo while this is modified.
2144 QScopedPointer<PathObject> original(new PathObject(*this));
2145
2146 // original_indices provides a mapping from this object's indices to the
2147 // equivalent original object indices in order to extract the relevant
2148 // parts of the original object later for cost calculation.
2149 // Note: curve handle indices may become incorrect, we don't need them.
2150 std::vector<MapCoordVector::size_type> original_indices;
2151 original_indices.resize(coords.size());
2152 for (std::size_t i = 0, end = coords.size(); i < end; ++i)
2153 original_indices[i] = i;
2154
2155 // A high value indicating an cost of deletion which is unknown.
2156 auto const undetermined_cost = std::numeric_limits<double>::max();
2157
2158 // This vector provides the costs of deleting individual nodes.
2159 std::vector<double> costs;
2160 costs.resize(coords.size(), undetermined_cost);
2161
2162 // The empty LineSymbol will not generate any renderables,
2163 // thus reducing the cost of update().
2164 // The temp object will be reused (but not reallocated) many times.
2165 LineSymbol empty_symbol;
2166 PathObject temp { &empty_symbol };
2167 temp.coords.reserve(10); // enough for two bezier edges.
2168
2169 for (auto part = path_parts.rbegin(); part != path_parts.rend(); ++part)
2170 {
2171 // Don't simplify parts which would get deleted.
2172 MapCoordVector::size_type minimum_part_size = part->isClosed() ? 4 : 3;
2173 auto minimumPartSizeReached = [&part, minimum_part_size]() -> bool
2174 {
2175 return (part->size() <= 7 && part->countRegularNodes() < minimum_part_size);
2176 };
2177
2178 auto minimum_cost_step = threshold / qMax(std::log2(part->size() / 64.0), 4.0);
2179 auto minimum_cost = 0.0;
2180
2181 // Delete nodes in max n runs (where n = max. part.end_index - part.start_index)
2182 for (auto run = part->last_index; run > part->first_index; --run)
2183 {
2184 if (minimumPartSizeReached())
2185 break;
2186
2187 // Determine the costs of node deletion
2188 {
2189 auto extract_start = part->first_index;
2190 auto index = part->nextCoordIndex(extract_start);
2191 while (index < part->last_index)
2192 {
2193 auto extract_end = part->nextCoordIndex(index);
2194 if (costs[index] == undetermined_cost)
2195 {
2196 temp.assignCoordinates(*this, extract_start, extract_end);
2197 temp.deleteCoordinate(index - extract_start, true, Settings::DeleteBezierPoint_RetainExistingShape);
2198 // Debug check: start and end coords of the extracts should be at the same position
2199 Q_ASSERT(coords[extract_start].isPositionEqualTo(temp.coords.front()));
2200 Q_ASSERT(coords[extract_end].isPositionEqualTo(temp.coords.back()));
2201 costs[index] = original->calcMaximumDistanceTo(original_indices[extract_start],
2202 original_indices[extract_end],
2203 &temp,
2204 0,
2205 temp.coords.size()-1);
2206 }
2207 extract_start = index;
2208 index = extract_end;
2209 }
2210
2211 if (costs[part->first_index] == undetermined_cost && extract_start < part->last_index)
2212 {
2213 if (part->isClosed())
2214 {
2215 auto extract_end = part->nextCoordIndex(part->first_index);
2216 temp.assignCoordinates(*this, extract_start, extract_end);
2217 temp.deleteCoordinate(part->last_index - extract_start, true, Settings::DeleteBezierPoint_RetainExistingShape);
2218 // Debug check: start and end coords of the extracts should be at the same position
2219 Q_ASSERT(coords[extract_start].isPositionEqualTo(temp.coords.front()));
2220 Q_ASSERT(coords[extract_end].isPositionEqualTo(temp.coords.back()));
2221 costs[part->first_index] = original->calcMaximumDistanceTo(original_indices[extract_start],
2222 original_indices[extract_end],
2223 &temp,
2224 0,
2225 temp.coords.size()-1);
2226 }
2227 else
2228 {
2229 costs[part->first_index] = threshold * 2;
2230 }
2231 }
2232
2233 costs[part->last_index] = undetermined_cost;
2234 }
2235
2236 // Find an upper bound for the acceptable cost in the current run
2237 auto cost_bound = threshold * 2;
2238 for (auto index = part->first_index; index < part->last_index; ++index)
2239 {
2240 if (costs[index] <= cost_bound)
2241 {
2242 cost_bound = costs[index];
2243 if (cost_bound <= minimum_cost)
2244 {
2245 // short-cut on minimal costs
2246 cost_bound = minimum_cost;
2247 break;
2248 }
2249 }
2250 }
2251
2252 if (cost_bound > threshold)
2253 break;
2254
2255 while (minimum_cost <= cost_bound)
2256 minimum_cost += minimum_cost_step;
2257 if (minimum_cost > threshold)
2258 minimum_cost = threshold;
2259
2260 // Start at the end, in order to shorten the vector quickly
2261 auto index = part->last_index;
2262 do
2263 {
2264 if (costs[index] <= cost_bound)
2265 {
2266 auto extract_start = part->prevCoordIndex(index);
2267 auto extract_end = part->nextCoordIndex(index);
2268 Q_ASSERT(original->coords[original_indices[extract_start]].isPositionEqualTo(coords[extract_start]));
2269 Q_ASSERT(original->coords[original_indices[extract_end]].isPositionEqualTo(coords[extract_end]));
2270 deleteCoordinate(index, true, Settings::DeleteBezierPoint_RetainExistingShape);
2271
2272 if (index > part->first_index)
2273 {
2274 // Deleted inner point
2275 auto new_extract_end = part->nextCoordIndex(extract_start);
2276 auto original_start = begin(original_indices) + extract_start;
2277 original_indices.erase(original_start + 1,
2278 original_start + 1 + (extract_end - new_extract_end));
2279 Q_ASSERT(original_indices.size() == coords.size());
2280 Q_ASSERT(original->coords[original_indices[extract_start]].isPositionEqualTo(coords[extract_start]));
2281 Q_ASSERT(original->coords[original_indices[new_extract_end]].isPositionEqualTo(coords[new_extract_end]));
2282
2283 costs.erase(begin(costs) + extract_start + 1,
2284 begin(costs) + extract_start + 1 + (extract_end - new_extract_end));
2285 Q_ASSERT(costs.size() == coords.size());
2286
2287 costs[extract_start] = undetermined_cost;
2288
2289 if (new_extract_end == part->last_index)
2290 {
2291 costs[part->first_index] = undetermined_cost;
2292 }
2293 else
2294 {
2295 costs[new_extract_end] = undetermined_cost;
2296 }
2297 Q_ASSERT(costs[part->last_index] > threshold);
2298
2299 index = part->prevCoordIndex(extract_start);
2300
2301 if (minimumPartSizeReached())
2302 break;
2303 }
2304 else if (part->isClosed())
2305 {
2306 Q_ASSERT(index == part->first_index);
2307
2308 // Deleted start point of closed path
2309 // This must match PathObject::deleteCoordinate
2310 extract_start = part->prevCoordIndex(part->last_index);
2311 auto original_begin = begin(original_indices);
2312 original_indices.erase(original_begin + part->first_index,
2313 original_begin + extract_end);
2314 original_indices.resize(coords.size());
2315 original_indices[part->last_index] = original_indices[part->first_index];
2316 Q_ASSERT(original->coords[original_indices[extract_start]].isPositionEqualTo(coords[extract_start]));
2317 Q_ASSERT(original->coords[original_indices[part->last_index]].isPositionEqualTo(coords[part->last_index]));
2318 Q_ASSERT(original->coords[original_indices[part->first_index]].isPositionEqualTo(coords[part->first_index]));
2319
2320 costs.erase(begin(costs) + part->first_index,
2321 begin(costs) + extract_end);
2322 costs.resize(coords.size(), undetermined_cost);
2323 costs[extract_start] = undetermined_cost;
2324 costs[part->first_index] = undetermined_cost;
2325 Q_ASSERT(costs[part->last_index] > threshold);
2326 }
2327 else
2328 {
2329 Q_ASSERT(index == part->first_index);
2330 }
2331 }
2332 else if (index != part->first_index)
2333 {
2334 --index;
2335 }
2336 }
2337 while (index != part->first_index);
2338 }
2339 }
2340
2341 bool removed_a_point = (coords.size() != original->coords.size());
2342 if (removed_a_point && undo_duplicate)
2343 {
2344 *undo_duplicate = original.take();
2345 }
2346
2347 return removed_a_point;
2348 }
2349
isPointOnPath(MapCoordF coord,qreal tolerance,bool treat_areas_as_paths,bool extended_selection) const2350 int PathObject::isPointOnPath(MapCoordF coord, qreal tolerance, bool treat_areas_as_paths, bool extended_selection) const
2351 {
2352 auto side_tolerance = tolerance;
2353 if (extended_selection && map && (symbol->getType() == Symbol::Line || symbol->getType() == Symbol::Combined))
2354 {
2355 // TODO: precalculate largest line extent for all symbols to move it out of this time critical method?
2356 side_tolerance = qMax(side_tolerance, symbol->calculateLargestLineExtent());
2357 }
2358
2359 auto contained_types = symbol->getContainedTypes();
2360 if ((contained_types & Symbol::Line || treat_areas_as_paths) && tolerance > 0)
2361 {
2362 update();
2363 for (const auto& part : path_parts)
2364 {
2365 const auto& path_coords = part.path_coords;
2366 auto size = path_coords.size();
2367 for (PathCoordVector::size_type i = 0; i < size - 1; ++i)
2368 {
2369 Q_ASSERT(path_coords[i].index < coords.size());
2370 if (coords[path_coords[i].index].isHolePoint())
2371 continue;
2372
2373 MapCoordF to_coord = coord - path_coords[i].pos;
2374 MapCoordF to_next = path_coords[i+1].pos - path_coords[i].pos;
2375 MapCoordF tangent = to_next;
2376 tangent.normalize();
2377
2378 auto dist_along_line = MapCoordF::dotProduct(to_coord, tangent);
2379 if (dist_along_line < -tolerance)
2380 continue;
2381
2382 if (dist_along_line < 0 && to_coord.lengthSquared() <= tolerance*tolerance)
2383 return Symbol::Line;
2384
2385 auto line_length = qreal(path_coords[i+1].clen) - qreal(path_coords[i].clen);
2386 if (line_length < 1e-7)
2387 continue;
2388
2389 if (dist_along_line > line_length + tolerance)
2390 continue;
2391
2392 if (dist_along_line > line_length && coord.distanceSquaredTo(path_coords[i+1].pos) <= tolerance*tolerance)
2393 return Symbol::Line;
2394
2395 auto right = tangent.perpRight();
2396
2397 auto dist_from_line = qAbs(MapCoordF::dotProduct(right, to_coord));
2398 if (dist_from_line <= side_tolerance)
2399 return Symbol::Line;
2400 }
2401 }
2402 }
2403
2404 // Check for area selection
2405 if ((contained_types & Symbol::Area) && !treat_areas_as_paths)
2406 {
2407 if (isPointInsideArea(coord))
2408 return Symbol::Area;
2409 }
2410
2411 return Symbol::NoSymbol;
2412 }
2413
isPointInsideArea(const MapCoordF & coord) const2414 bool PathObject::isPointInsideArea(const MapCoordF& coord) const
2415 {
2416 update();
2417 bool inside = false;
2418 for (const auto& part : path_parts)
2419 {
2420 if (part.isPointInside(coord))
2421 inside = !inside;
2422 }
2423 return inside;
2424 }
2425
calcMaximumDistanceTo(MapCoordVector::size_type start_index,MapCoordVector::size_type end_index,const PathObject * other,MapCoordVector::size_type other_start_index,MapCoordVector::size_type other_end_index) const2426 double PathObject::calcMaximumDistanceTo(
2427 MapCoordVector::size_type start_index,
2428 MapCoordVector::size_type end_index,
2429 const PathObject* other,
2430 MapCoordVector::size_type other_start_index,
2431 MapCoordVector::size_type other_end_index) const
2432 {
2433 update();
2434
2435 Q_ASSERT(other_start_index == 0);
2436 Q_ASSERT(other_end_index == other->coords.size()-1);
2437
2438 if (end_index < start_index)
2439 {
2440 const auto part = findPartForIndex(start_index);
2441 Q_ASSERT(part->isClosed());
2442 Q_ASSERT(end_index >= part->first_index);
2443 Q_ASSERT(end_index <= part->last_index);
2444 auto d1 = calcMaximumDistanceTo(start_index, part->last_index, other, other_start_index, other_end_index);
2445 auto d2 = calcMaximumDistanceTo(part->first_index, end_index, other, other_start_index, other_end_index);
2446 return qMax(d1, d2);
2447 }
2448
2449 const float test_points_per_mm = 2;
2450
2451 float max_distance_sq = 0.0;
2452 for (const auto& part : path_parts)
2453 {
2454 if (part.first_index <= end_index && part.last_index >= start_index )
2455 {
2456 const auto& path_coords = part.path_coords;
2457
2458 auto pc_start = std::lower_bound(begin(path_coords), end(path_coords), start_index, PathCoord::indexLessThanValue);
2459 auto pc_end = std::lower_bound(pc_start, end(path_coords), end_index, PathCoord::indexLessThanValue);
2460 if (pc_end == end(path_coords))
2461 {
2462 if (pc_end != pc_start)
2463 --pc_end;
2464 }
2465
2466 PathCoord path_coord;
2467 float distance_sq = 0.0;
2468 other->calcClosestPointOnPath(pc_start->pos, distance_sq, path_coord, other_start_index, other_end_index);
2469 max_distance_sq = qMax(max_distance_sq, distance_sq);
2470
2471 for (auto pc = pc_start; pc != pc_end; ++pc)
2472 {
2473 auto next_pc = pc + 1;
2474 double len = next_pc->clen - pc->clen;
2475 MapCoordF direction = next_pc->pos - pc->pos;
2476
2477 int num_test_points = qMax(1, qRound(len * test_points_per_mm));
2478 for (int p = 1; p <= num_test_points; ++p)
2479 {
2480 MapCoordF point = pc->pos + direction * ((float)p / num_test_points);
2481 other->calcClosestPointOnPath(point, distance_sq, path_coord, other_start_index, other_end_index);
2482 max_distance_sq = qMax(max_distance_sq, distance_sq);
2483 }
2484 }
2485 }
2486 }
2487
2488 return sqrt(max_distance_sq);
2489 }
2490
makeIntersectionAt(double a,double b,const PathCoord & a0,const PathCoord & a1,const PathCoord & b0,const PathCoord & b1,PathPartVector::size_type part_index,PathPartVector::size_type other_part_index)2491 PathObject::Intersection PathObject::Intersection::makeIntersectionAt(double a, double b, const PathCoord& a0, const PathCoord& a1, const PathCoord& b0, const PathCoord& b1, PathPartVector::size_type part_index, PathPartVector::size_type other_part_index)
2492 {
2493 PathObject::Intersection new_intersection;
2494 new_intersection.coord = MapCoordF(a0.pos.x() + a * (a1.pos.x() - a0.pos.x()), a0.pos.y() + a * (a1.pos.y() - a0.pos.y()));
2495 new_intersection.part_index = part_index;
2496 new_intersection.length = a0.clen + a * (a1.clen - a0.clen);
2497 new_intersection.other_part_index = other_part_index;
2498 new_intersection.other_length = b0.clen + b * (b1.clen - b0.clen);
2499 return new_intersection;
2500 }
2501
normalize()2502 void PathObject::Intersections::normalize()
2503 {
2504 std::sort(begin(), end());
2505 erase(std::unique(begin(), end()), end());
2506 }
2507
calcAllIntersectionsWith(const PathObject * other,PathObject::Intersections & out) const2508 void PathObject::calcAllIntersectionsWith(const PathObject* other, PathObject::Intersections& out) const
2509 {
2510 update();
2511 other->update();
2512
2513 const double epsilon = 1e-10;
2514 const double zero_minus_epsilon = 0 - epsilon;
2515 const double one_plus_epsilon = 1 + epsilon;
2516
2517 for (size_t part_index = 0; part_index < path_parts.size(); ++part_index)
2518 {
2519 const PathPart& part = path_parts[part_index];
2520 auto path_coord_end_index = part.path_coords.size() - 1;
2521 for (auto i = PathCoordVector::size_type { 1 }; i <= path_coord_end_index; ++i)
2522 {
2523 // Get information about this path coord
2524 bool has_segment_before = (i > 1) || part.isClosed();
2525 MapCoordF ingoing_direction;
2526 if (has_segment_before && i == 1)
2527 {
2528 Q_ASSERT(path_coord_end_index >= 1);
2529 ingoing_direction = part.path_coords[path_coord_end_index].pos - part.path_coords[path_coord_end_index - 1].pos;
2530 ingoing_direction.normalize();
2531 }
2532 else if (has_segment_before)
2533 {
2534 Q_ASSERT(i >= 1 && i < part.path_coords.size());
2535 ingoing_direction = part.path_coords[i-1].pos - part.path_coords[i-2].pos;
2536 ingoing_direction.normalize();
2537 }
2538
2539 bool has_segment_after = (i < path_coord_end_index) || part.isClosed();
2540 MapCoordF outgoing_direction;
2541 if (has_segment_after && i == path_coord_end_index)
2542 {
2543 Q_ASSERT(part.path_coords.size() > 1);
2544 outgoing_direction = part.path_coords[1].pos - part.path_coords[0].pos;
2545 outgoing_direction.normalize();
2546 }
2547 else if (has_segment_after)
2548 {
2549 Q_ASSERT(i < path_coord_end_index);
2550 outgoing_direction = part.path_coords[i+1].pos - part.path_coords[i].pos;
2551 outgoing_direction.normalize();
2552 }
2553
2554 // Collision state with other object at current other path coord
2555 bool colliding = false;
2556 // Last known intersecting point.
2557 // This is valid as long as colliding == true and entered as intersection
2558 // when the next segment suddenly is not colliding anymore.
2559 Intersection last_intersection;
2560
2561 for (size_t other_part_index = 0; other_part_index < other->path_parts.size(); ++other_part_index)
2562 {
2563 const PathPart& other_part = other->path_parts[part_index]; /// \todo FIXME: part_index or other_part_index ???
2564 auto other_path_coord_end_index = other_part.path_coords.size() - 1;
2565 for (auto k = PathCoordVector::size_type { 1 }; k <= other_path_coord_end_index; ++k)
2566 {
2567 // Test the two line segments against each other.
2568 // Naming: segment in this path is a, segment in other path is b
2569 const PathCoord& a0 = part.path_coords[i-1];
2570 const PathCoord& a1 = part.path_coords[i];
2571 const PathCoord& b0 = other_part.path_coords[k-1];
2572 const PathCoord& b1 = other_part.path_coords[k];
2573 MapCoordF b_direction = b1.pos - b0.pos;
2574 b_direction.normalize();
2575
2576 bool first_other_segment = (k == 1);
2577 if (first_other_segment)
2578 {
2579 colliding = isPointOnSegment(a0.pos, a1.pos, b0.pos);
2580 if (colliding && !other_part.isClosed())
2581 {
2582 // Enter intersection at start of other segment
2583 bool ok;
2584 double a = parameterOfPointOnLine(a0.pos.x(), a0.pos.y(), a1.pos.x() - a0.pos.x(), a1.pos.y() - a0.pos.y(), b0.pos.x(), b0.pos.y(), ok);
2585 Q_ASSERT(ok);
2586 out.push_back(Intersection::makeIntersectionAt(a, 0, a0, a1, b0, b1, part_index, other_part_index));
2587 }
2588 }
2589
2590 bool last_other_segment = (k == other_path_coord_end_index);
2591 if (last_other_segment)
2592 {
2593 bool collision_at_end = isPointOnSegment(a0.pos, a1.pos, b1.pos);
2594 if (collision_at_end && !other_part.isClosed())
2595 {
2596 // Enter intersection at end of other segment
2597 bool ok;
2598 double a = parameterOfPointOnLine(a0.pos.x(), a0.pos.y(), a1.pos.x() - a0.pos.x(), a1.pos.y() - a0.pos.y(), b1.pos.x(), b1.pos.y(), ok);
2599 Q_ASSERT(ok);
2600 out.push_back(Intersection::makeIntersectionAt(a, 1, a0, a1, b0, b1, part_index, other_part_index));
2601 }
2602 }
2603
2604 double denominator = a0.pos.x()*b0.pos.y() - a0.pos.y()*b0.pos.x() - a0.pos.x()*b1.pos.y() - a1.pos.x()*b0.pos.y() + a0.pos.y()*b1.pos.x() + a1.pos.y()*b0.pos.x() + a1.pos.x()*b1.pos.y() - a1.pos.y()*b1.pos.x();
2605 if (denominator == 0)
2606 {
2607 // Parallel lines, calculate parameters for b's start and end points in a and b.
2608 // This also checks whether the lines are actually on the same level.
2609 bool ok;
2610 double b_start = 0;
2611 double a_start = parameterOfPointOnLine(a0.pos.x(), a0.pos.y(), a1.pos.x() - a0.pos.x(), a1.pos.y() - a0.pos.y(), b0.pos.x(), b0.pos.y(), ok);
2612 if (!ok)
2613 {
2614 if (colliding) out.push_back(last_intersection);
2615 colliding = false;
2616 continue;
2617 }
2618 double b_end = 1;
2619 double a_end = parameterOfPointOnLine(a0.pos.x(), a0.pos.y(), a1.pos.x() - a0.pos.x(), a1.pos.y() - a0.pos.y(), b1.pos.x(), b1.pos.y(), ok);
2620 if (!ok)
2621 {
2622 if (colliding) out.push_back(last_intersection);
2623 colliding = false;
2624 continue;
2625 }
2626
2627 // Cull ranges
2628 if (a_start < zero_minus_epsilon && a_end < zero_minus_epsilon)
2629 {
2630 if (colliding) out.push_back(last_intersection);
2631 colliding = false;
2632 continue;
2633 }
2634 if (a_start > one_plus_epsilon && a_end > one_plus_epsilon)
2635 {
2636 if (colliding) out.push_back(last_intersection);
2637 colliding = false;
2638 continue;
2639 }
2640
2641 // b overlaps somehow with a, check if we have to enter one or two collisions
2642 // (provided the incoming / outgoing tangents are not parallel!)
2643 if (!colliding)
2644 {
2645 if (a_start <= 0)
2646 {
2647 // b comes in over the start of a
2648 Q_ASSERT(a_end >= 0);
2649
2650 // Check for parallel tangent case
2651 if (!has_segment_before || MapCoordF::dotProduct(b_direction, ingoing_direction) < 1 - epsilon)
2652 {
2653 // Enter intersection at a=0
2654 double b = b_start + (0 - a_start) / (a_end - a_start) * (b_end - b_start);
2655 out.push_back(Intersection::makeIntersectionAt(0, b, a0, a1, b0, b1, part_index, other_part_index));
2656 }
2657
2658 colliding = true;
2659 }
2660 else if (a_start >= 1)
2661 {
2662 // b comes in over the end of a
2663 Q_ASSERT(a_end <= 1);
2664
2665 // Check for parallel tangent case
2666 if (!has_segment_after || -1 * MapCoordF::dotProduct(b_direction, outgoing_direction) < 1 - epsilon)
2667 {
2668 // Enter intersection at a=1
2669 double b = b_start + (1 - a_start) / (a_end - a_start) * (b_end - b_start);
2670 out.push_back(Intersection::makeIntersectionAt(1, b, a0, a1, b0, b1, part_index, other_part_index));
2671 }
2672
2673 colliding = true;
2674 }
2675 else
2676 Q_ASSERT(false);
2677 }
2678 if (colliding)
2679 {
2680 if (a_end > 1)
2681 {
2682 // b goes out over the end of a
2683 Q_ASSERT(a_start <= 1);
2684
2685 // Check for parallel tangent case
2686 if (!has_segment_after || MapCoordF::dotProduct(b_direction, outgoing_direction) < 1 - epsilon)
2687 {
2688 // Enter intersection at a=1
2689 double b = b_start + (1 - a_start) / (a_end - a_start) * (b_end - b_start);
2690 out.push_back(Intersection::makeIntersectionAt(1, b, a0, a1, b0, b1, part_index, other_part_index));
2691 }
2692
2693 colliding = false;
2694 }
2695 else if (a_end < 0)
2696 {
2697 // b goes out over the start of a
2698 Q_ASSERT(a_start >= 1);
2699
2700 // Check for parallel tangent case
2701 if (!has_segment_before || -1 * MapCoordF::dotProduct(b_direction, ingoing_direction) < 1 - epsilon)
2702 {
2703 // Enter intersection at a=0
2704 double b = b_start + (0 - a_start) / (a_end - a_start) * (b_end - b_start);
2705 out.push_back(Intersection::makeIntersectionAt(1, b, a0, a1, b0, b1, part_index, other_part_index));
2706 }
2707
2708 colliding = false;
2709 }
2710 else
2711 {
2712 // b stops in the middle of a.
2713 // Remember last known colliding point, the endpoint of b.
2714 last_intersection = Intersection::makeIntersectionAt(a_end, 1, a0, a1, b0, b1, part_index, other_part_index);
2715 }
2716 }
2717
2718 // Check if there is a collision at the endpoint
2719 Q_ASSERT(colliding == (a_end >= 0 && a_end <= 1));
2720 }
2721 else
2722 {
2723 // Non-parallel lines, calculate intersection parameters and check if in range
2724 double a = +(a0.pos.x()*b0.pos.y() - a0.pos.y()*b0.pos.x() - a0.pos.x()*b1.pos.y() + a0.pos.y()*b1.pos.x() + b0.pos.x()*b1.pos.y() - b1.pos.x()*b0.pos.y()) / denominator;
2725 if (a < zero_minus_epsilon || a > one_plus_epsilon)
2726 {
2727 if (colliding) out.push_back(last_intersection);
2728 colliding = false;
2729 continue;
2730 }
2731
2732 double b = -(a0.pos.x()*a1.pos.y() - a1.pos.x()*a0.pos.y() - a0.pos.x()*b0.pos.y() + a0.pos.y()*b0.pos.x() + a1.pos.x()*b0.pos.y() - a1.pos.y()*b0.pos.x()) / denominator;
2733 if (b < zero_minus_epsilon || b > one_plus_epsilon)
2734 {
2735 if (colliding) out.push_back(last_intersection);
2736 colliding = false;
2737 continue;
2738 }
2739
2740 // Special case for overlapping (cloned / traced) polylines: check if b is parallel to adjacent direction.
2741 // If so, set colliding to true / false without entering an intersection because the other line
2742 // simply continues along / comes from the path of both polylines instead of intersecting.
2743 double dot = -1;
2744 if (has_segment_before && a <= 0 + epsilon)
2745 {
2746 // Ingoing direction
2747 dot = MapCoordF::dotProduct(b_direction, ingoing_direction);
2748 if (b <= 0 + epsilon)
2749 dot = -1 * dot;
2750 }
2751 else if (has_segment_after && a >= 1 - epsilon)
2752 {
2753 // Outgoing direction
2754 dot = MapCoordF::dotProduct(b_direction, outgoing_direction);
2755 if (b >= 1 - epsilon)
2756 dot = -1 * dot;
2757 }
2758 if (dot >= 1 - epsilon)
2759 {
2760 colliding = (b > 0.5);
2761 continue;
2762 }
2763
2764 // Enter the intersection
2765 last_intersection = Intersection::makeIntersectionAt(a, b, a0, a1, b0, b1, part_index, other_part_index);
2766 out.push_back(last_intersection);
2767 colliding = (b == 1);
2768 }
2769 }
2770 }
2771 }
2772 }
2773 }
2774
setCoordinate(MapCoordVector::size_type pos,const MapCoord & c)2775 void PathObject::setCoordinate(MapCoordVector::size_type pos, const MapCoord& c)
2776 {
2777 Q_ASSERT(pos < getCoordinateCount());
2778
2779 const PathPart& part = *findPartForIndex(pos);
2780 if (part.isClosed() && pos == part.last_index)
2781 pos = part.first_index;
2782 coords[pos] = c;
2783 if (part.isClosed() && pos == part.first_index)
2784 setClosingPoint(part.last_index, c);
2785
2786 setOutputDirty();
2787 }
2788
addCoordinate(MapCoordVector::size_type pos,const MapCoord & c)2789 void PathObject::addCoordinate(MapCoordVector::size_type pos, const MapCoord& c)
2790 {
2791 Q_ASSERT(pos <= coords.size());
2792
2793 if (coords.empty())
2794 {
2795 coords.emplace_back(c);
2796
2797 path_parts.clear();
2798 path_parts.emplace_back(*this, 0, 0);
2799 }
2800 else
2801 {
2802 auto part = findPartForIndex(qMin(pos, coords.size() - 1));
2803 coords.insert(coords.begin() + pos, c);
2804 partSizeChanged(part, +1);
2805
2806 if (pos == part->first_index)
2807 {
2808 if (part->isClosed())
2809 setClosingPoint(part->last_index, c);
2810 }
2811 else if (pos == part->last_index)
2812 {
2813 Q_ASSERT(pos > part->first_index);
2814 coords[pos-1].setClosePoint(false);
2815 coords[pos-1].setHolePoint(false);
2816 }
2817 }
2818 setOutputDirty();
2819 }
2820
addCoordinate(const MapCoord & c,bool start_new_part)2821 void PathObject::addCoordinate(const MapCoord& c, bool start_new_part)
2822 {
2823 if (coords.empty())
2824 {
2825 addCoordinate(0, c);
2826 }
2827 else if (!start_new_part)
2828 {
2829 addCoordinate(coords.size(), c);
2830 }
2831 else
2832 {
2833 auto index = coords.size();
2834 coords.push_back(c);
2835 path_parts.emplace_back(*this, index, index);
2836 setOutputDirty();
2837 }
2838 }
2839
deleteCoordinate(MapCoordVector::size_type pos,bool adjust_other_coords,int delete_bezier_point_action)2840 void PathObject::deleteCoordinate(MapCoordVector::size_type pos, bool adjust_other_coords, int delete_bezier_point_action)
2841 {
2842 const auto part = findPartForIndex(pos);
2843 const auto coords_begin = begin(coords);
2844
2845 auto num_coords = part->size();
2846 if (num_coords <= 7)
2847 {
2848 auto num_regular_points = part->countRegularNodes();
2849
2850 if (num_regular_points <= 2 && (pos == part->first_index || pos == part->last_index))
2851 {
2852 // Too small, must delete
2853 deletePart(findPartIndexForIndex(pos));
2854 return;
2855 }
2856
2857 if (num_regular_points == 3 && num_coords == 4)
2858 {
2859 // A closed path of three straight segments, to be opened
2860 Q_ASSERT(part->isClosed());
2861 if (pos == part->last_index)
2862 pos = part->first_index;
2863
2864 // Move pos to end_index-1, than erase last two coords
2865 std::rotate(coords_begin + part->first_index, coords_begin + pos + 1, coords_begin + part->last_index);
2866 coords.erase(coords_begin + part->last_index - 1, coords_begin + part->last_index+1);
2867 partSizeChanged(part, -2);
2868 coords[part->last_index].setHolePoint(true);
2869 return;
2870 }
2871 }
2872
2873 auto prev_index = part->prevCoordIndex(pos);
2874 if (prev_index < pos && coords[prev_index].isCurveStart() && prev_index + 3 > pos)
2875 {
2876 // Delete curve handles
2877 coords[prev_index].setCurveStart(false);
2878 coords.erase(coords_begin + prev_index + 1, coords_begin + prev_index + 3);
2879 partSizeChanged(part, -2);
2880 return;
2881 }
2882
2883 if (pos == part->first_index || pos == part->last_index)
2884 {
2885 if (part->isClosed())
2886 {
2887 // Make start/end point an inner point before removing
2888 auto middle_offset = coords[part->first_index].isCurveStart() ? 3 : 1;
2889 std::rotate(coords_begin + part->first_index, coords_begin + part->first_index + middle_offset, coords_begin + part->last_index);
2890 setClosingPoint(part->last_index, coords[part->first_index]);
2891 pos = part->last_index - middle_offset;
2892 prev_index = part->prevCoordIndex(pos);
2893 }
2894 else
2895 {
2896 auto start = std::min<VirtualPath::size_type>(pos, prev_index + 1);
2897 auto end = std::max<VirtualPath::size_type>(pos + 1, part->nextCoordIndex(pos));
2898 coords.erase(coords_begin + start, coords_begin + end);
2899 partSizeChanged(part, start - end);
2900 coords[part->last_index].setHolePoint(true);
2901 coords[part->last_index].setCurveStart(false);
2902 return;
2903 }
2904 }
2905
2906 // Delete inner point
2907 Q_ASSERT(pos > part->first_index);
2908 Q_ASSERT(pos < part->last_index);
2909 if (!(coords[prev_index].isCurveStart() && coords[pos].isCurveStart()))
2910 {
2911 // Bezier curve on single side
2912 if (coords[pos].isCurveStart())
2913 coords[prev_index].setCurveStart(true);
2914
2915 coords.erase(coords_begin + pos);
2916 partSizeChanged(part, -1);
2917 }
2918 else
2919 {
2920 // Bezier curves on both sides
2921 if (adjust_other_coords)
2922 {
2923 // Adjust handle positions to preserve the original shape somewhat
2924 // (of course, in general it is impossible to do this)
2925 prepareDeleteBezierPoint(pos, delete_bezier_point_action);
2926 }
2927
2928 coords.erase(begin(coords) + pos - 1, begin(coords) + pos + 2);
2929 partSizeChanged(part, -3);
2930 }
2931 }
2932
prepareDeleteBezierPoint(MapCoordVector::size_type pos,int delete_bezier_point_action)2933 void PathObject::prepareDeleteBezierPoint(MapCoordVector::size_type pos, int delete_bezier_point_action)
2934 {
2935 const MapCoord& p0 = coords[pos - 3];
2936 MapCoord& p1 = coords[pos - 2];
2937 const MapCoord& p2 = coords[pos - 1];
2938 const MapCoord& q0 = coords[pos];
2939 const MapCoord& q1 = coords[pos + 1];
2940 MapCoord& q2 = coords[pos + 2];
2941 const MapCoord& q3 = coords[pos + 3];
2942
2943 double pfactor, qfactor;
2944 if (delete_bezier_point_action == Settings::DeleteBezierPoint_ResetHandles)
2945 {
2946 double target_length = BEZIER_HANDLE_DISTANCE * p0.distanceTo(q3);
2947 pfactor = target_length / qMax(p0.distanceTo(p1), 0.01);
2948 qfactor = target_length / qMax(q3.distanceTo(q2), 0.01);
2949 }
2950 else if (delete_bezier_point_action == Settings::DeleteBezierPoint_RetainExistingShape)
2951 {
2952 calcBezierPointDeletionRetainingShapeFactors(p0, p1, p2, q0, q1, q2, q3, pfactor, qfactor);
2953
2954 if (qIsInf(pfactor))
2955 pfactor = 1;
2956 if (qIsInf(qfactor))
2957 qfactor = 1;
2958
2959 calcBezierPointDeletionRetainingShapeOptimization(p0, p1, p2, q0, q1, q2, q3, pfactor, qfactor);
2960
2961 double minimum_length = 0.01 * p0.distanceTo(q3);
2962 pfactor = qMax(minimum_length / qMax(p0.distanceTo(p1), 0.01), pfactor);
2963 qfactor = qMax(minimum_length / qMax(q3.distanceTo(q2), 0.01), qfactor);
2964 }
2965 else
2966 {
2967 Q_ASSERT(delete_bezier_point_action == Settings::DeleteBezierPoint_KeepHandles);
2968 pfactor = 1;
2969 qfactor = 1;
2970 }
2971
2972 MapCoordF p0p1 = MapCoordF(p1) - MapCoordF(p0);
2973 p1 = MapCoord(p0.x() + pfactor * p0p1.x(), p0.y() + pfactor * p0p1.y());
2974
2975 MapCoordF q3q2 = MapCoordF(q2) - MapCoordF(q3);
2976 q2 = MapCoord(q3.x() + qfactor * q3q2.x(), q3.y() + qfactor * q3q2.y());
2977 }
2978
clearCoordinates()2979 void PathObject::clearCoordinates()
2980 {
2981 coords.clear();
2982 path_parts.clear();
2983 setOutputDirty();
2984 }
2985
assignCoordinates(const PathObject & proto,MapCoordVector::size_type first,MapCoordVector::size_type last)2986 void PathObject::assignCoordinates(const PathObject& proto, MapCoordVector::size_type first, MapCoordVector::size_type last)
2987 {
2988 Q_ASSERT(last < proto.coords.size());
2989
2990 auto part = proto.findPartForIndex(first);
2991 Q_ASSERT(part == proto.findPartForIndex(last));
2992
2993 coords.clear();
2994 if (last >= first)
2995 coords.reserve(last - first + 1);
2996 else
2997 coords.reserve(last - part->first_index + part->last_index - first + 2);
2998
2999 for (auto i = first; i != last; )
3000 {
3001 MapCoord new_coord = proto.coords[i];
3002 new_coord.setClosePoint(false);
3003 new_coord.setHolePoint(false);
3004 coords.push_back(new_coord);
3005
3006 ++i;
3007 if (i > part->last_index)
3008 {
3009 i = part->first_index;
3010 if (part->isClosed() && i != last)
3011 coords.erase(coords.end()-1);
3012 }
3013
3014 }
3015 // Add last point
3016 MapCoord new_coord = proto.coords[last];
3017 new_coord.setCurveStart(false);
3018 new_coord.setClosePoint(false);
3019 new_coord.setHolePoint(true);
3020 coords.push_back(new_coord);
3021
3022 path_parts.clear();
3023 path_parts.emplace_back(*this, 0, coords.size()-1);
3024
3025 setOutputDirty();
3026 }
3027
updatePathCoords() const3028 void PathObject::updatePathCoords() const
3029 {
3030 auto part_start = VirtualPath::size_type { 0 };
3031 for (auto& part : path_parts)
3032 {
3033 part.first_index = part_start;
3034 part.last_index = part.path_coords.update(part_start);
3035 part_start = part.last_index+1;
3036 }
3037 }
3038
recalculateParts()3039 void PathObject::recalculateParts()
3040 {
3041 setOutputDirty();
3042
3043 path_parts.clear();
3044 if (!coords.empty())
3045 {
3046 MapCoordVector::size_type start_index = 0;
3047 auto last_index = coords.size()-1;
3048
3049 for (MapCoordVector::size_type i = 0; i <= last_index; ++i)
3050 {
3051 if (coords[i].isHolePoint())
3052 {
3053 path_parts.emplace_back(*this, start_index, i);
3054 start_index = i+1;
3055 }
3056 else if (coords[i].isCurveStart())
3057 {
3058 i += 2;
3059 }
3060 }
3061
3062 if (start_index <= last_index)
3063 path_parts.emplace_back(*this, start_index, last_index);
3064 }
3065 }
3066
setClosingPoint(MapCoordVector::size_type index,const MapCoord & coord)3067 void PathObject::setClosingPoint(MapCoordVector::size_type index, const MapCoord& coord)
3068 {
3069 auto& out_coord = coords[index];
3070 out_coord = coord;
3071 out_coord.setCurveStart(false);
3072 out_coord.setHolePoint(true);
3073 out_coord.setClosePoint(true);
3074 }
3075
updateEvent() const3076 void PathObject::updateEvent() const
3077 {
3078 updatePathCoords();
3079 }
3080
createRenderables(ObjectRenderables & output,Symbol::RenderableOptions options) const3081 void PathObject::createRenderables(ObjectRenderables& output, Symbol::RenderableOptions options) const
3082 {
3083 symbol->createRenderables(this, path_parts, output, options);
3084 }
3085
3086
3087 // ### PointObject ###
3088
PointObject(const Symbol * symbol)3089 PointObject::PointObject(const Symbol* symbol)
3090 : Object(Object::Point, symbol, { MapCoord{0,0} })
3091 {
3092 Q_ASSERT(!symbol || (symbol->getType() == Symbol::Point));
3093 }
3094
3095 PointObject::PointObject(const PointObject& /*proto*/) = default;
3096
3097
duplicate() const3098 PointObject* PointObject::duplicate() const
3099 {
3100 return new PointObject(*this);
3101 }
3102
copyFrom(const Object & other)3103 void PointObject::copyFrom(const Object& other)
3104 {
3105 if (&other == this)
3106 return;
3107
3108 Object::copyFrom(other);
3109 const PointObject* point_other = other.asPoint();
3110 if (getSymbol() && getSymbol()->isRotatable())
3111 setRotation(point_other->getRotation());
3112 }
3113
3114
setPosition(qint32 x,qint32 y)3115 void PointObject::setPosition(qint32 x, qint32 y)
3116 {
3117 coords[0].setNativeX(x);
3118 coords[0].setNativeY(y);
3119 setOutputDirty();
3120 }
3121
setPosition(const MapCoord & coord)3122 void PointObject::setPosition(const MapCoord& coord)
3123 {
3124 coords[0] = coord;
3125 }
3126
setPosition(const MapCoordF & coord)3127 void PointObject::setPosition(const MapCoordF& coord)
3128 {
3129 coords[0].setX(coord.x());
3130 coords[0].setY(coord.y());
3131 setOutputDirty();
3132 }
3133
getCoordF() const3134 MapCoordF PointObject::getCoordF() const
3135 {
3136 return MapCoordF(coords.front());
3137 }
3138
getCoord() const3139 MapCoord PointObject::getCoord() const
3140 {
3141 return coords.front();
3142 }
3143
transform(const QTransform & t)3144 void PointObject::transform(const QTransform& t)
3145 {
3146 if (t.isIdentity())
3147 return;
3148
3149 auto& coord = coords.front();
3150 const auto p = t.map(MapCoordF{coord});
3151 coord.setX(p.x());
3152 coord.setY(p.y());
3153 setOutputDirty();
3154 }
3155
3156
intersectsBox(const QRectF & box) const3157 bool PointObject::intersectsBox(const QRectF& box) const
3158 {
3159 return box.contains(QPointF(coords.front()));
3160 }
3161
3162
3163 } // namespace OpenOrienteering
3164